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Abstract 

This  thesis  considers  the  problem  of  positioning  a  very  long,  vertical  tow  cable  in  the  ocean. 
Drill  strings  and  tethers  that  support  remotely-operated  vehicles  are  examples  of  systems 
for  which  position  control  of  the  lower  endpoint  is  important.  Consistent  with  ctirrent 
practice,  we  assume  that  the  cable  can  be  maneuvered  only  by  vessel  motions  at  the  top, 
and  that  it  is  tracked  only  at  the  endpoints. 

A  derivation  of  the  equations  of  cable  motion  is  given,  followed  by  an  analysis  of  the 
nonlinear  frequency  response  of  the  plant,  using  the  method  of  harmonic  balances  and  a  per¬ 
turbation  approach.  The  perturbation  results  are  extended  to  give  a  closed-form  estimate 
of  the  closed-loop  limit-cycling  behavior.  The  basic  control  approach  is  to  consider  input 
preshaping  and  regulation  designs  separately,  merging  them  to  form  tracking  controllers. 
The  primary  preshaping  part  is  a  frequency-domain  dynamic  inversion  based  on  recent  re¬ 
sults  in  robotics;  it  requires  the  construction  of  a  finite-dimensional  model  of  the  plant, 
and  works  for  in-plane  and  coupled  out-of-plane  motions,  as  well  as  some  other  distributed- 
pareimeter  physical  systems.  For  regulator  design,  a  number  of  observer-based  approaches 
are  considered,  ranging  from  standard  linear  loopshaping  to  an  approximately  optimal  non¬ 
linear  control  law  with  nonlinear  observation.  Our  solution  to  the  nonlinear  optimal  control 
problem  is  novel  in  the  sense  that  it  can  accommodate  plant  models  of  arbitrary  order,  a 
necessity  for  distributed  plants. 

The  preshaping  techniques  are  verified  with  full-scale  results  from  experiments  in  the 
ocean  with  2000  meters  of  cable,  and  the  various  closed-loop  methods  are  compared  based 
on  scsde-model  laboratory  tests. 

Thesis  Supervisor:  Dana  R.  Yoerger 

Title:  Associate  Scientist,  Woods  Hole  Oceanographic  Institution 
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Chapter  1 


Introduction 


1.1  Motivation 

Towed  onderwater  cables  have  been  an  important  component  of  ocean  systems  for  many 
decades.  In  early  oceanographic  work  the  cable  was  a  simple  steel  wire  for  lowering  sampling 
instruments,  corers,  dredging  tools,  and,  less  commonly,  bathyspheres  to  specified  depths. 
More  recently,  these  cables  have  been  constructed  of  state-of-the-axt  materials,  and  are  often 
capable  of  supporting  very  high  communication  bandwidths  through  fiber  optic  technology. 
The  cables  have  also  grown  in  diameter  to  enable  the  transmission  of  substantial  electrical 
and  mechanical  power.  With  these  advances,  cables  rem^  an  extremely  effective  way  to 
reach  the  deep  seafloor. 

High-tension  cables  form  an  important  subclass  for  towed  systems.  The  tension  is  due  to 
the  weight  of  the  mass  at  the  bottom  (referred  to  in  the  rest  of  the  thesis  as  “the  towfish”), 
and  to  the  weight  of  the  cable  itself.  High-tension  cables  are  popular  because  the  weight 
keeps  them  essentially  vertical  under  the  surface  vessel.  There  is  little  danger  of  hitting  the 
cable  with  the  props,  and  approximate  positioning  of  the  towfish  is  accomplished  merely 
by  positioning  the  surface  vessel.  Currents  acting  in  the  water  column  cause  the  cable  to 
billow  only  a  small  amount. 

There  are  many  situations  in  which  the  towfish  must  be  positioned  with  an  accuracy  of 
ten  meters  or  less.  For  a  deep  system,  one  solution  is  to  put  thrusters  on  the  towfish,  with  the 
idea  that  the  towfish  can  reach  all  points  in  a  circle  of  given  radius  while  the  surface  vessel  is 
stationary.  The  size  of  this  “footprint”  is  naturally  governed  by  the  thrust  capability  of  the 
towfish,  and  by  the  water  weight  of  the  system.  Another  approach  to  precise  positioning, 


13 


the  one  which  is  developed  in  this  thesis,  is  to  incoiporate  our  knowledge  of  the  cable’s 
dynamic  response,  and  position  the  towfish  solely  by  maneuvering  the  surface  vessel.  This 
eliminates  the  need  to  exert  large  forces  at  the  towfish  and  frees  the  towfish  from  the  power 
and  space  requirements  of  thrusters.  For  very  heavy  systems,  the  reqtiired  thruster  forces 
could  be  exceedingly  large. 

Figure  1-1  shows  three  high-tension  ocean  towing  applications  where  accurate  position¬ 
ing  is  important  and  towfish  thrusters  are  not  used.  In  a)  is  drawn  a  deep-sea  drilling 
deployment.  The  object  is  to  position  the  drill  bit  assembly  at  the  lower  end  accurately 
enough  so  that  an  existing  hole  can  be  reentered.  As  drawn,  a  cone  of  several  meters’  radiiis 
is  often  used  to  align  the  riser  with  the  hole.  The  riser  may  be  forty  centimeters  or  more  in 
diameter  and  weigh  more  than  one  thousand  Newtons  per  meter  submerged.  The  surface 
vessel  is  maneuvered  with  dynamic  positioning,  and  the  bit  assembly  is  tracked  using  a 
short-baseline  acoustic  net. 

In  part  b)  of  Figure  1-1,  a  two- vehicle  ROV  system  for  exploration,  survey,  or  inspection 
is  shown.  The  vertical  cable  is  under  high  tension,  while  the  tether  coimecting  the  towfish 
to  the  ROV  is  of  very  low  tension.  The  presence  of  the  heavy  towfish  and  the  low-tension 
cable  isolates  the  ROV  from  heave  motions  of  the  surface  vessel,  making  the  ROV  free  to 
move  within  the  tether  radius  using  its  own  thrusters.  In  practice,  it  is  the  towfish  that  must 
be  positioned  within  the  tether  radius  of  the  ROV,  since  the  ROV  is  interacting  with  the 
seafloor.  The  ARGO/JASON  system,  developed  at  the  Deep  Submergence  Laboratory  at 
the  Woods  Hole  Oceanographic  Institution,  is  similar  to  that  drawn.  The  ARGO/JASON 
vertical  cable  is  steel  and  1.7  centimeters  in  diameter,  and  the  towfish  weighs  approximately 
4500  Newtons  in  water.  The  towfish  is  typically  tracked  using  a  long-baseline  acoustic  net. 

An  envisioned  scenario  for  ocean  mining  is  drawn  in  Figure  1-lc).  Although  ocean 
miuiTig  has  not  been  pursued  on  a  large  scale  yet,  there  is  great  interest,  since  large  quantities 
of  ferromanganese  nodules  have  been  located  on  the  seafloor  [3].  The  proposed  systems 
would  be  capable  of  recovering  some  5000  tons  of  ore  per  day.  One  recent  design  calls  for 
a  large  bottom-crawling  collector  which  grinds  the  ore  and  separates  it  from  the  gangue;  a 
pumping  unit,  drawn  here  as  the  towfish,  provides  the  means  to  move  the  material  along  a 
full-depth  vertical  pipeline.  This  mining  system  is  a  cross  between  the  systems  of  Figures  1- 
la)  and  1-lb)  in  that  the  main  cable  is  extremely  heavy  and  of  a  large  diameter,  and  the 
motions  of  the  pumping  unit  are  slaved  to  those  of  the  collector. 
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Figiue  1-1:  Three  high-tension  ocean  towing  applications:  a)  Borehole  reentry  with  a  drill 
string,  b)  ROV  system,  c)  Ocean  mining. 
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Accurate  positioning  of  the  towhsh  is  required  in  all  of  these  applications.  This  the¬ 
sis  develops  appropriate  techniques  in  open-loop  and  closed-loop  senses,  without  reqtiiring 
thrusters  at  the  towfish. 


1.2  Research  Objectives 

The  plant  under  consideration  is  a  difficult  one  to  control,  when  tackled  with  the  standard 
design  procedures.  The  system  is  governed  by  partial  differential  equations,  and  can  have 
significant  geometric  nonlinearities  in  addition  to  quadratic  drag.  Furthermore,  we  can  only 
reasonably  measure  the  position  of  a  single  or  a  few  points  along  the  member,  and  exert 
control  on  the  system  only  at  the  top  end. 

It  is  the  goal  of  this  thesis  to  investigate  a  number  of  approaches,  based  on  classical 
control  theory,  that  can  be  applied  to  the  positioning  problem.  We  note  at  the  outset  that 
very  few  formal  guarantees  can  be  made,  because  an  extreme  simplification  of  the  plant  is 
used.  Nonetheless,  it  it  einticipated  that  designs  based  on  the  simple  models  will  yield  good 
properties  from  an  operational  viewpoint.  The  objectives  of  this  research  are  as  follows: 

1.  Provide  an  analytical  approximation  of  the  lateral  frequency  response  of  towed  vehicle 
systems.  Such  an  analysis  will  give  some  quantification  of  the  expected  behavior 
through  simple  calculations. 

2.  Develop  a  technique  for  improved  dynamic  response  during  standard  point-to-point 
maneuvers,  as  woiild  be  required  for  an  efficient  change  of  position  for  ARGO/JASON. 
This  procedure  will  improve  upon  the  very  slow  settling  times  experienced  during 
conventional  maneuvers. 

3.  Develop  a  technique  for  computing  vessel  trajectories  that  lead  to  precision  path- 
following  for  the  lower  endpoint  of  the  member,  as  would  be  required  for  a  sonar 
survey  following  a  specific  search  pattern.  Given  the  very  long  time-delay  and  the 
nonlinear  nature  of  the  cable,  we  hope  to  **sharpen  up”  the  general  response. 

4.  Develop  feedback  control  algorithms  for  station-keeping  and  path-following  in  the 
presence  of  unmodeled  dynamics,  sensor  noise,  and  physical  disturbances  such  as 
currents.  Feedback  control  is  crucial  in  drill-hole  reentry,  for  example,  where  a  given 
watch  circle  of  several  meters  must  be  maintained  with  the  drill  bit. 
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5.  Verify  the  above  procedures  experimentally.  Since  no  model  is  absolutely  correct, 
experimental  verification  lends  more  credibility  to  the  theoretical  developments.  In 
addition,  experiments  axe  likdy  to  reveal  the  shortcomings  of  the  design:  ^f^rbiiiques, 
and  to  point  to  directions  of  further  work. 

1.3  Outline  of  the  Thesis 

Chapter  2  begins  with  a  brief  but  complete  derivation  of  the  governing  equations  of  motion 
for  a  cable-towfish  system.  The  general  nonlinear  frequency  response  is  approximated  using 
harmonic  balances  and  a  perturbation  technique;  the  effects  of  ambient  current  are  included 
in  the  analysis.  The  perturbation  formulation  is  extended  to  provide  a  closed-form  estimate 
of  the  limit-cycling  behavior  of  a  cable  system  under  linear  closed-loop  control.  Comparisons 
with  simulation  data  are  presented  for  all  the  approximations,  as  verification. 

The  open-loop  positioning  problem  is  treated  in  Chapter  3.  The  rather  ad  hoc  method 
of  overshooting  for  point-to-point  moves  is  briefly  discussed,  and  then  a  ftiU  derivation  of 
a  dynamic  inversion  scheme  is  given.  This  scheme  has  a  proof  of  convergence,  and  may  be 
applied  to  a  range  of  distributed-parameter  problems.  Full-scale  experiments  are  described, 
which  verify  the  point-to-point  and  trajectory-foUowing  techniques. 

Chapter  4  discusses  a  number  of  important  issues  relating  to  regulator  design,  which  is 
to  be  observer-based.  The  basic  properties  of  the  standard  estimators  and  control  laws  are 
listed,  and  a  new  approach  for  approximating  the  solution  to  the  optimal  nonlinear  feedback 
problem  is  derived.  This  approximate  solution  is  novel  and  practical  in  the  sense  that  it  is 
nondynamicai,  does  not  require  a  look-up  table,  and  can  be  applied  to  models  of  arbitrary 
order.  Furthermore,  the  tools  of  its  derivation  may  be  useful  in  controUer  design  for  other 
plants  with  high  order  and  nonlinearities  with  low-order  polynomial  approximations.  At 
the  close  of  the  chapter,  a  set  of  proposed  compensators  is  specified. 

The  application  of  closed-loop  regulation  and  tracking  is  treated  in  Chapter  5.  Some 
of  the  practical  issues  that  arise  in  implementation  of  the  theory  are  discussed,  such  as  the 
explicit  calculation  of  tolerable  modeling  errors  in  the  different  feedback  strategies.  This  is 
a  particularly  interesting  problem  for  nonlinear  feedback  because,  as  will  be  seen,  the  imcer- 
tainty  must  be  expressed  at  the  plant  input.  Two  scale-model  experiments  are  described  and 
data  are  provided  which  demonstrate  the  relative  merits  of  each  control  approach.  These 
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data  cover  both  regulation  and  tracking;  the  latter  is  accomplished  by  simply  layering  the 
regulators  of  Chapter  4  on  top  of  the  precalculated  trajectories  of  Chapter  3. 

Chapter  6  gives  conclusions  and  recommendations,  including  an  assessment  of  the  prac¬ 
tical  advantages  and  disadvantages  of  the  various  estimation  and  control  methods. 

1.4  Related  Literature 

The  dynamic  responses  of  tethers  and  cables  have  been  modeled  accurately  for  a  long 
time.  References  with  relevance  to  ocean  applications  include  Sanders  [87],  Ablow  and 
Schecter  [2],  Triantafyllou  [97],  Bliek  [19],  Burgess  [23],  and  Howell  [56].  The  dynamics  of 
underwater  vehicles  are  also  well  known;  the  seminal  formulas  have  been  given  by  Abkowitz 
[1],  while  the  standard  set  of  hydrodynamic  terms  has  been  provided  by  Feldman  [42]. 
The  response  of  towed  underwater  cable  systems  to  vessel  motions  has  been  discussed  by 
Mudie  and  Iverson  [71],  and  Chapman  [27].  These  papers  are  primarily  concerned  with  the 
transmission  of  vessel  sway  motions  down  the  cable,  and  with  the  cable  response  during 
turns.  An  experimental  study  of  half-turn  maneuvers  was  carried  out  by  Guerch  [52], 
verifying  the  model  developed  in  the  paper  by  Iverson  and  Mudie  [57]. 

The  idea  of  preshaping  actuator  trajectories  for  complex  motions  of  noncolocated  sys¬ 
tems  has  been  used  extensively  in  robotics  applications.  One  of  the  originators  was  Smith 
[93],  whose  “posicast”  approach  utilized  direct  knowledge  of  the  system’s  frequency  re¬ 
sponse.  The  important  papers  of  Bayo  and  his  coworkers  [9],  [11],  [12],  [88]  cover  a  range  of 
approaches  in  robotics;  the  idea  is  to  specify  the  desired  trajectory  of  the  endpoint  and  then 
explicitly  compute  the  required  motor  torque  to  achieve  it.  Singer,  Seering,  and  Singhose 
[89]  [90]  took  a  similar  approach,  with  the  goal  of  eliminating  residual  vibrations  in  a  robotic 
arm  at  the  end  of  a  move.  In  contrast  to  Bayo’s  work,  these  papers  use  simpler  plant  models 
with  only  several  natural  modes. 

Although  experienced  operators  at  sea  have  long  been  able  to  execute  remarkably  ef¬ 
fective  trajectories,  the  idea  of  actually  precomputing  vessel  trajectories  to  effect  specific 
towfish  motions  is  only  fifteen  years  old.  An  early  work  is  Pa\il  and  Soler  [80],  where  com¬ 
plete  trajectories  for  fast  point-to-point,  in-plane  maneuvers  were  generated,  using  iteration 
and  finite-element  modeling.  In  a  different  approach,  predictive  dynamic  modeling  was  used 
in  real  time  by  Mudie  and  Hastens  [72]  to  specify  vessel  trajectories  that  bring  the  towfish 
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closest  (in  the  horizontal  plane)  to  a  specified  target.  These  researchers  conducted  sea  tests 
with  3300  meters  of  cable,  and  reported  an  average  miss  distance  of  65  meters  coming  out 
of  a  turn  of  approximately  two-kilometer  radius.  Cheng  [29]  has  further  investigated  the 
use  of  predictive  models;  his  approach  provides  a  full  display  for  the  operator,  who  is  to 
make  the  decisions  on  course  changes.  The  technique  was  verified  using  simulations  and 
human  subjects. 

As  with  the  open-loop  problem,  feedback  control  of  distributed  parameter  systems  has 
been  explored  primarily  in  the  area  of  robotics.  There  are  a  very  large  number  of  papers 
on  this  topic,  so  we  list  only  a  few  of  the  most  important  works  here.  Modal  control  for 
generic  systems  was  discussed  by  Balas  [6],  and  Cannon  and  Schmitz  [25]  provided  one  of 
the  first  application  papers,  using  the  Linear  Quadratic  Gaussian/Loop  Transfer  Recovery 
(LQG/LTR)  method  on  a  finite-element  model  of  a  single-link  arm.  Sakawa  et  al.  [86] 
addressed  the  flexible  robot  feedback  problem  with  a  modal  approach,  the  sliding  mode  was 
used  by  Yeung  and  Chen  [104],  and  an  adaptive  approach  with  prediction  was  developed 
by  Cetinkunt  and  Wu  [26]. 

Vincet  et  al.  [101]  recognized  that  to  reposition  noncolocated  linear  mass-spring  sys¬ 
tems,  it  often  makes  sense  to  decompose  the  control  law  into  two  parts.  The  first  part 
is  a  saturated-input  high-gain  control  law  to  execute  the  gross  movement,  followed  by  a 
passive  control  law  to  eliminate  residual  vibrations.  In  contrast  to  the  previous  references 
on  feedback  control,  full-state  measurement  or  estimation  is  not  required  here.  Lyapunov 
arguments  (e.g.,  [92])  indeed  verify  that  any  passive  control  law  applied  at  the  boundary  of 
such  a  system  (which  could  be  nonlinear)  leads  to  asymptotic  stability.  These  observations 
are  consistent  with  those  of  the  paper  by  Liu  et  al.  [62],  who  showed  that  a  system  of  vibrat¬ 
ing  strings  can  be  exponentially  stabilized  by  placing  linear  dampers  at  each  or  both  ends. 
This  idea  has  been  applied  to  vibration  suppression  in  robots  and  structures  by  Paden  et 
al.  [78]  (passive  joint  controllers),  and  Hagedom  [53]  (active  absorption  of  traveling  waves). 

In  this  thesis,  the  idea  of  splitting  the  control  law  into  two  parts  will  be  used.  The 
feedforward  part  is  related  to  Bayo’s  work,  and  the  feedback  laws  are  model-based,  in  a 
robust  control  setting.  For  the  linear  single-input /single-output  (SISO)  feedback  problem, 
one  would  typicaUy  use  classical  loopshaping  ideas  for  this  (e.g.,  [38]),  while  for  m\ilti- 
input/multi-ouput  systems  (MIMO),  the  LQG/LTR  approach  [39],  [94]  is  popular.  The 
synthesis  of  nonlinear  estimators  and  nonlinear  regulators  was  discussed  by  Gnmberg  and 
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Athans  [50],  [51],  who  provided  a  loop  recovery  procedure  much  LTR,  called  Loop  Operator 
Recovery  (LOR).  This  result  will  be  used  heavily  in  Chapter  4,  and  actually  implemented 
in  Chapter  5.  ^ 

The  nonlinear  estimation  problem  has  been  recently  discussed  in  the  survey  paper  by  Mi- 
sawa  and  Hedrick  [68].  Among  the  most  popular  of  techniques,  the  constant  gain  extended 
Kalman  filter  [84]  and  the  full  extended  Kalman  filter  (e.g.,  [44])  have  good  properties  that 
will  be  discussed  in  Chapter  4.  On  the  control  side,  although  the  LQR  has  some  intrinsic 
robustness  when  applied  to  nonlinear  systems  [82],  a  nonlinear  optimal  control  approxima¬ 
tion  wUl  be  worked  out  as  well.  The  statement  of  the  optimal  control  problem  for  nonlinear 
systems  can  be  foimd  in  the  book  by  Athans  and  Falb  [5],  and  the  properties  of  optimal 
controllers  with  quadratic  cost  were  disctissed  by  Moylan  and  Anderson  [70],  Glad  [46],  and 
Tsitsiklis  and  Athans  [99].  Actual  solutions  of  the  Hamilton-Jacobi-Bellman  (HJB)  equa¬ 
tion  have  become  possible  with  today’s  faster  computers;  methods  of  computation  have 
been  discussed  by  Beaman  [13],  Dwyer  [40],  O’Sullivan  and  Sain  [77],  Luus  [66],  [67],  and 
Nagurka  and  Yen  [73].  The  analytical  approach  of  this  thesis  is  based  on  the  work  of  Yoshida 
and  Loparo  [105],  and  is  similar  to  that  of  Willemstein  [103].  The  latter  approach  has  been 
applied  to  satellite  control  by  Dwyer  and  Sen  [41]  and  Dabbous  and  Ahmed  [32]. 

The  paper  by  Triantafyllou  and  Grosenbaugh  [98]  is  one  of  the  very  few  applications  of 
feedback  control  to  the  towing  problem.  The  authors  adopted  a  linear  second-order  lumped 
mass  model  of  the  cable  system,  with  a  pure  time  delay,  and  applied  the  Smith  regulator 
idea  [93]  to  the  LQG/LTR  methodology.  A  high-order  model  developed  by  Burgess  and 
Triantafyllou  [24]  was  used  for  verification,  in  an  in-plane  example.  The  application  of  LQG 
control  to  ocean  mining  problems  has  been  discussed  by  Okano  [76],  who  reported  that  in 
large  pipe  deployments,  the  frequencies  of  vortex-induced  vibration  [64]  are  usually  well- 
separated  from  the  first  five  natural  lateral  modes  of  vibration  for  the  pipe.  Finally,  the 
effect  of  atmospheric  drag  on  slewing  control  of  flexible  space  robots  has  been  investigated 
by  Juang  and  Horta  [58].  Written  from  the  viewpoint  of  extrapolating  the  response  in  a 
vacuum  from  experiments  in  air,  these  researchers  found  that  the  primary  effect  of  the  drag 
was  to  eliminate  high-frequency  vibrational  modes.  We  expect  a  similar  trend  to  be  true  in 
ocean  cables,  and  this  point  will  be  investigated  further  in  Section  2.6.2. 


‘It  is  expected  that  the  more  recent  approaches  such  as  Rao  (e.g.,  [43],  [47])  and  the  p-synthesis  approach 
(e.g.,  [37])  might  lead  to  similar  results,  although  as  yet  there  are  no  tractable  nonlinear  extensions. 
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Chapter  2 


Cable  Dynamic  Equations  and 
Analysis 


2.1  Introduction 

In  order  to  address  our  control  aims  in  the  context  of  underwater  towed  systems,  it  is  impor¬ 
tant  to  model  the  physical  system  as  accurately  as  possible,  while  maintaining  a  formulation 
that  is  reasonable  for  controller  design  purposes.  This  chapter  provides  a  derivation  of  the 
relevant  equations  of  motion  for  cables,  and  discusses  the  approximations  involved  in  mak¬ 
ing  a  model  of  the  plant  for  control.  Because  of  quadratic  drag,  even  the  simplified  models 
have  a  complex  dynamics,  so  we  give  two  approximate  analyses  of  the  lateral  frequency 
response  of  a  vertical  imderwater  cable.  The  frequency-response  results  are  then  extended 
to  predict  the  limit-cycle  behavior  of  a  closed-loop  system,  with  a  linear  control  law  and  a 
single  measurement.  Limit  cycles  are  commonly  observed  in  nonlinear  control  systems  of 
this  type. 


2.2  Cable  Equations 

We  begin  with  the  full  three-dimensional  coupled  equations,  and  show  how  they  may  be 
simplified.  The  equations  without  bending  stiffness  have  been  derived  in  the  lagrangian 
coordinate  system  by  Bliek  [19]  for  example,  and  bending  stiffness  was  subsequently  treated 
in  Howell  [56].  The  following  development  is  based  on  these  sources.  An  initial  assumption 
is  that  the  cable  is  inextensible,  consistent  with  the  fact  that  axial  vibrations  usually  have 
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a  small  effect  on  the  lateral  response  of  a  vertical  cable.  In  more  general  cases,  where  there 
is  significant  curvature  in  the  cable,  the  geometric  coupling  can  be  important. 

2.2.1  Reference  Systems 

The  problem  may  be  set  in  a  lagrangian  frame  with  the  tangential,  normeJ,  and  binormal 
unit  vectors  written  as  [t,  n,  &].  The  Eiiler  angle  convention  provides  a  way  to  relate  quan¬ 
tities  in  the  inertial  frame  [X,  Y,  Z]  to  those  in  the  local  frame.  The  ordering  of  angles  is 
specified  as  shown  in  Figure  2-1; 

1.  Rotate  by  <{>  around  Z;  X  —*  X'  and  Y  Y'. 

2.  Rotate  by  6  around  F';  X'  — ►  X"  and  Z  —*  Z\  This  brings  X”  in  line  with  £. 

3.  Rotate  by  xp  around  X"  to  bring  F'  — *  n  and  Z'  -♦  b. 

Substantial  derivatives  (/?(■))  in  moving  reference  frames  involve  the  Euler  angles;  let 
C3  be  the  rotation  vector  (the  time  rate  of  change  of  the  cable  element  orientation),  and  0 
be  the  Darboux  vector  (the  spatijd  rate  of  change  of  the  cable  element  orientation).  Also, 
let  t  denote  time,  and  s  be  the  lagrangian  cable  coordinate.  Then  for  an  arbitrary  vector 
G  :=  Gii  -t-  G2n  -|-  Gab  in  the  local  frame,  we  have: 


22 


(2.1) 


DG  ac  _  = 


where  u>  =  wii  +  wjn  +  u^b,  ft  =  fti<  +  ftjh  +  flab,  and 


di{j  d4>  .  ^ 

— - —stnB^ 

dt  dt 

dB  d4> 

—cosxjf  +  —cos6sinil>y 
dt  dt 

d<f>  dB 

—cos$costj>  —  —sinrp, 

dt  dt 


dfj)  d<j)  .  n 

dJ  - 

dO  ,  d<f>  ^  , 
—cosw  +  -r-coawtnv’, 
ds  ds 

d<t>  .  ,  dB  .  , 

—cosBcosyf - —stnrp. 

ds  ds 


2.2.2  Force  and  Moment  Equations 

From  Figure  2-2,  Newton’s  Law  gives: 


+  max9  =  ^  +  ilxF  +  i:R,  (2.5) 

dt  os 

where  m  is  the  cable  mass  per  unit  length,  V  is  the  local  velocity  vector,  F  is  the  internal 
force  vector,  and  'LR  is  the  the  aggregate  external  force  vector.  Normally,  LA  consists  of 
gravity,  drag,  and  added  mass  terms. 

To  allow  for  bending  stiffness  terms,  we  must  also  write  out  the  moment  equations.  Let 
M  be  the  local  internal  moment  vector;  if  the  cross-section  is  circular,  solid,  and  homoge¬ 
neous,  the  moments  can  be  written  as  follows: 
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ds 


M  = 


GJ  0  Q 
Q  El  Q 
Q  Q  El 


fi. 


(2.6) 


Here,  G  is  the  shear  modudus,  J  is  the  polar  moment  of  inertia,  E  is  Young’s  modulus,  and 
I  is  the  bending  moment  of  inertia.  We  define  a  matrix 


H  = 


J  0  0 
0/0 
0  0/ 


so  that  the  moment  equation  is 


da  dM  _  (It 

H—  +  <5  X  {HCS)  =  i^  +  nxM  +  (?  +  df-Li2+  —  XT. 
ot  09  ds 


(2.7) 


Here,  f  is  the  position  vector  of  the  cable  element,  and  Q  is  the  applied  moment  to  the 
element. 


2.2.3  Simplifications 

We  neglect  distributed  externa]  moments  Q,  and  note  that  in  the  limit  as  df  goes  to  zero, 
drx  EE  goes  to  zero  also,  and  ^  goes  to  t.  Thus,  we  have  a  simpler  form  for  Equation  2.7: 


24 


(2.8) 


dd  .  dM  ^  ^  ^ 

+  w  X  (Hd)  =  ^  +  (lxM  +  txT. 
dt  '  '  d» 

Howell  [56]  carried  out  a  nondimensional  analysis  which  showed  that  certain  terms  of 
the  force  and  moment  vector  equations  can  be  neglected  based  on  the  physics  of  the  cable 
system.  These  points  will  be  taken  into  account  here;  they  are  outlined  below. 

1.  The  ratio  of  inertial  moment  to  bending  moment  in  rotation  is  on  the  order  of  or 
10~^X  for  steel  cables,  and  lower  for  synthetic  cables.  Thus,  we  may  usually  neglect 
the  rotational  inertia  terms.  Here,  pc  is  the  density  of  the  cable,  g  is  the  acceleration 
due  to  gravity,  and  L  is  the  length  of  the  cable. 

2.  With  no  rotational  inertia  terms,  the  torsion  on  the  length  must  be  constant.  If  the 
cable  is  torque*balanced  and  there  are  no  applied  moments  at  the  boundaries,  then 
the  torsion  effects  are  zero.  We  may  thus  simplify  matters  by  setting  the  Euler  angle 
V-  =  0. 

3.  In  nondimensionalized  variables,  the  shear  forces  F2  and  F3  are  )>  '*^^**'e 

Wo  is  the  weight  in  water  per  unit  length  of  the  cable,  and  W  is  the  in-water  weight 
of  the  towfish.  The  importance  of  bending  stiffness  terms  can  be  judged  on  this 
quantity;  high  tension  and  large  length  generally  make  the  bending  moments  small, 
and  possibly  negligible.  For  the  current  work,  we  expect  that  this  ratio  will  be  small, 
so  the  bending  stiffness  terms  are  included  through  the  derivations  only.  An  exception 
is  that,  in  Section  5.2.2,  which  deals  with  model  scaling  of  deep  systems,  this  ratio  is 
calculated  to  verify  that  the  scaling  is  adequate  with  respect  to  bending  stiffness. 


2.2.4  Compatibility  Relations 

In  addition  to  force  (Equation  2.5)  and  moment  relations  (Equation  2.8),  a  cable  is  also 
subject  to  certain  kinematic  constraints,  referred  to  as  the  compatibility  relations.  We  have 


D_£f_ 

DtDs  ~  DsDt' 


and  since  ^  =  V  and  ^  =  f,  this  is  just 


U>  Xt  =  -r—  -1-  n  X  V. 

as 


(2.9) 


(2.10) 
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2.2.5  Governing  Equations 

The  full  set  of  governing  equations  consist  of  three  force  equations  (2.5),  two  bending 
equations  (2.8  with  no  torsional  dynamics),  and  three  compatibility  relations  (2.10).  Let 
[u,  V,  w]  be  velocities  in  the  tangential,  normal,  and  binormal  directions,  respectively.  The 
bending  equations  are  not  prognostic,  so  the  tmlcnown  variables  are  u,  v,  w,  Fi,  and  0. 


2.2.6  Applied  Forces 

We  now  consider  the  applied  forces  in  the  lagrangian  coordinate  system.  They  consist  of 
gravity,  drag,  and  added  mass:  EF  =  EF(u,  w, (j>,9)  =  Rg-\-  Rd-\-  Ram-  The  separate  parts 
are  as  follows. 

1.  Rg  =  [-Woain<l>cos9,  —WoCos4>,  —WoSin<f>sin9],  where  Y  is  taken  to  be  upwards  in  the 
inertial  frame. 

2.  S.J  =  [-^pdCtV.  I  u  I,  -^pdCdV-s/v^  +  w*,  -  \pdCjw-\/v^  +  tr*],  where  p  is  the  den¬ 
sity  of  water,  Cj  and  Ct  are  the  normal  and  tangential  drag  coefficients  of  the  cable, 
respectively,  and  d  is  the  diameter  of  the  cable.  We  assume  that  the  same  drag  coef¬ 
ficient  is  valid  in  both  the  normal  and  binormal  directions,  using  Morison’s  approach 
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181]. 

3.  Ram  =  [0,  -nia^,  where  is  the  added  mass  per  unit  length.  Added  mass 

in  the  axi2d  direction  is  neglected. 

In  situations  with  currents,  the  velocities  in  Rj  above  should  be  taken  as  the  velocities 
relative  to  the  ambient  flow  field.  In  addition,  we  have  neglected  hydrodynamic  lift  forces. 

2.2.7  In-Plane  Motions 

The  above  equations  collapse  into  a  set  of  four  equations  if  only  in-plane  motions  are 
considered.  There  is  only  one  Euler  angle  in  this  case;  the  unknowns  are  u,v,<f>,  and  Fi. 
For  convenience,  from  here  on  the  tension  is  denoted  by  a  T  instead  of  Fi .  We  have 


du  d4> 

m— - Tnv-;r- 

dt  dt 

.  .dv  b4> 

dv 

ds 

du 

ds 


dT^„.d<i>dH  •  a  1  I  1 
=  T~  ~  ~w^s4>-~f>dCjv  \vl, 

d£_  ^ 

dt  “as’ 


(2.12) 


Equations  2.12  represent  the  complete  equations  of  motion  for  an  ineztensible  in-plane 
cable.  They  are  augmented  by  appropriate  boundary  conditions,  three  at  each  end.  For 
OUT  applications,  the  top  end  and  the  bottom  end  are  pinned  (no  transmitted  moment), 
and  the  top  and  bottom  motions  are  given  by  the  vessel  forcing  and  the  towfish  dynamics, 
respectively. 


2.3  General  Behavior  of  Vertical  Cables 

Vertical  cables  in  water  are  governed  by  a  damped  wave  equation,  which  may  or  may  not 
take  into  account  bending  stiflhess.  This  fact  indicates  that  a  transfer  function  which  takes 
motions  of  the  top  end  to  motions  of  the  bottom  end  has  a  pxire  delay.  If  there  is  large 
bending  stiffness,  then  the  system  can  be  non-minimtun  phase  as  well.  These  two  properties 
are  notorious  for  making  controller  design  difficult.  However,  simulations  show  that  the 
extreme  forces  of  quadratic  drag  make  it  impossible  to  transmit  any  but  the  smallest  of 
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motions  all  the  way  down  an  underwater  cable.  The  cable  thus  serves  as  a  very  strong  filter 
of  motions,  whose  gain  increases  as  the  amplitude  of  motion  decreases.  This  situation  is 
illustrated  in  Figure  2-3,  where  simulation  of  Equation  2.38  (see  Section  2.5.1)  has  been 
used  to  create  an  amplitude-dependent  Bode  plot.  It  is  clear  that  when  the  excitation  is 
very  small,  the  gain  and  phase  approach  those  of  an  undamped  wave  equation  transfer 
function.  As  the  excitation  is  increased,  the  gain  peaks  are  reduced  and  shifted  to  lower 
frequencies,  and  the  phase  is  smoothed  and  reduced  (toward  — oo)  on  the  frequency  axis. 
These  observations  are  consistent  with  the  well-known  behavior  of  linear  systems  when 
damping  is  increased.  As  the  system  is  in  fact  nonlinear,  the  next  two  sections  discuss  how 
this  dependence  of  dynamic  response  upon  the  amplitude  of  excitation  may  be  characterized 
further.  Understanding  the  nonlinear  behavior  of  this  system  is  important  from  a  control 
point  of  view;  for  example,  the  attenuation  which  accompanies  large  excitations  is  stabilizing 
by  virtue  of  the  small-gain  theorem  [33],  but  the  loss  of  phase  that  comes  with  these  large 
motions  can  destabilize  linear  control  systems. 

Before  treating  the  approximations,  we  note  that  Figure  2-3  above,  and  all  of  the  other 
calculations  for  the  chapter,  are  based  on  the  nominal  system  described  in  Table  2.3. 

2.4  Harmonic  Balance 

The  first  approximation  technique  we  consider  is  a  variation  of  the  well-known  harmonic 
balance  or  equivalent  linearization  procedure  to  find  the  in-plane  dynamics  of  mooring 
systems  [48].  This  approach  allows  one  to  include  the  effects  of  static  curvature  (caused  by 
currents  perhaps),  and  elongation  of  the  cable  under  dynamic  loading.  As  such,  and  for 
completeness,  the  derivations  which  follow  include  these  effects,  even  though  Equations  2.12 
do  not. 


2.4.1  Cable  Equations 

We  consider  the  expansion  T  T  +  T,  where  f  is  the  static  tension,  and  T  is  the  dynamic 
tension.  Similarly,  we  let  ^  ~  ^  -f  The  static  position  of  the  cable  is  given  uniquely  by 
so  we  let  p  and  q  denote  the  axial  and  lateral  dynamic  deflections,  respectively.  The  cable 
coordinate  system  is  oriented  in  the  inertial  frame  as  shown  in  Fig\ire  2-4. 

To  compute  the  static  configuration  in  the  presence  of  a  current  U  (positive  in  the  X- 
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log(fre4)uency,  nd/tec) 


log(fre<)uency,  tad/sec) 


ise  of  the  nominal  system:  gain  (top)  and  phase 
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cable  length 

cable  diameter 

cable  water  weight /length 

towfish  water  weight 

cable  mass/length 

cable  added  mass  (lateral)/length 

effective  horizontal  mass  of  towfish 

effective  vertical  mass  of  towfish 

effective  rotational  inertia  of  towfish 

water  density 

cable  normal  drag  coefficient 
cable  tangential  drag  coefficient 
drag  coefficient  of  towfish 
area  of  towfish  vertical  face 
area  of  towfish  horizontal  face 
linear  rotational  damping  of  towfish 
distance  between  center  of  gravity  and 
cable  connection  point 
cable  elasticity 


L 

1000  m 

d 

0.02  meters 

Wo 

15  Nfm 

W 

4000  AT 

m 

1.6  kgfm 

rria 

0.4  kgfm 

Mp  +  MaP 

500  kg 

Mq  +  MaQ 

500  kg 

J 

100  kg  •  m* 

p 

1025  kg/m^ 

Cd 

2 

Ct 

0.01 

2 

Ap 

1  m* 

Aq 

2  m* 

Bi 

100  N  •  m/{m/s) 

Ac 

1  m 

E 

200  X  10»  JV/m* 

Table  2.1:  Nominal  system  parameters. 


direction),  one  can  use  a  drag  model  of  arbitrary  complexity,  as  the  static  shape  has  to  be 
computed  only  once.  For  example,  a  popular  form  that  can  induce  lift  is 

£),  (2.13) 

This  particular  form  can  lead  to  decreased  static  tensions  during  towing  because  the  lift 
forces  carry  some  of  the  weight,  while  the  drag  forces  become  largely  axial  and  therefore  are 
controlled  by  Ct-  As  written,  these  static  drag  terms  are  influenced  by  elongation,  in  that 
stretching  causes  a  decrease  in  diameter  through  Poisson’s  ratio.  In  addition,  elongation  will 
always  show  up  in  the  representation  of  the  static  cable  configuration  in  world  coordinates. 

In  the  dynamics,  we  assume  that  the  quadratic  drag  acting  upon  the  cable  can  be 
treated  as  a  linear  damping  with  spatial  dependence.  In  particular,  the  describing  function 
of  ((s)sin(ijt)  I  sm(uit)  |  is  ((s)^sin(u>t);  the  cases  of  nonzero  c\irrents  will  be  described 


=  ^pdCfUcos^  I  Ucos^  I  ^1  + 
Rd,2  =  —^pdCdUsin^  \  Usin^  \  fl 
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in  Section  2.4.3.  Under  the  assumptions  of  small  dynamic  angle  and  precomputed  static 
configuration,  the  d}mamic  equations  of  Section  2.2.7  are: 


df 

ds 


ds 

ds 


d^p  -  -  1  _  dp  I  dp  I 

^pdCjU\U\-w^sin$^-f^, 

f  ^  d^ 

EA'^^ds' 


(2.14) 

(2.15) 


(2.16) 

(2.17) 


where  we  have  written  the  compatibility  relations  in  terms  of  positions  instead  of  velocities. 
Adopting  linear  damping  coefficients  bp(s),  bg(s)  which  are  to  be  foimd  iteratively,  one  can 
write  this  in  the  frequency  domain  in  matrix  form: 


WoCos^  — mu>*  +  tu;&p(s) 
« 

0  0 


0 

f 

y^(-(m  +  ma)u;2  +  tw6,(s)) 

P 

0 

.  ?  . 

or  =  A(s)y.  The  constant  term  of  Equation  2.15  transforms  to  a  delta  function,  so  it 
has  no  bearing  on  the  frequency  response.  We  let  the  cable  be  discretized  into  n  nodes, 
numbered  1  at  the  towfish  connection  point,  and  n  at  the  surface  point,  where  there  is  an 
imposed  motion.  The  unstretched  distance  between  nodes  is  ds.  Thus,  we  may  approximate 
the  evolution  of  y  on  the  length  as 


yi+i  -  Vi  ^  ^(^+i)^+i  +  i4(si)yi 
ds  ~  2 


(2.18) 


W+i 


/-yA(si+,)]  '[/+yA(50]!7i:=5,y. 


Then  y„  is  related  to  j/i  as  follows: 


(2.19) 
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Vn  =  Bn~lBn-2  *  '‘Biffi  :=  0yi. 


(2.20) 


This  represents  four  equations  in  six  unknowns:  Ti,  T„,  and  pi  and  91. 

2.4.2  Towfish  Dynamics 

The  required  additional  information  is  found  by  looking  at  the  towfish  dynamics;  there  is  a 
relationship  between  the  tension  and  angle  at  the  connection  point,  and  the  motions  of  the 
towfish.  With  the  variables  P,  Q,  and  i  defined  as  shown  in  Figure  2*4,  the  basic  equations 
for  the  towfish  are  as  follows: 


d^P 

(Mp  + 

(Mq  + 

^  dt^  dt 


1  dP 

Ticoa{<f>i  -  #)  -  Wsini  -  ■^pCd,towfi$hAp— 

1  dQ 

Txsin{4>\  -i)-  Wcoai  -  -pCd^towfithAq-^ 


dP 
dt  ’ 

aQ 

dt  ’ 


-T\AeCOs{<f>x  —  #). 


(2.21) 

(2.22) 

(2.23) 


The  static  solution  requires  that  -  #  =  y,  so  that  Equation  2.22  can  be  solved  alone  to 
obtain  #.  Again,  one  could  use  a  sophisticated  drag  law  for  the  vehicle  in  the  steady  state; 
for  example,  the  static  drag  force  in  the  JT-direction  can  be  approximated  as 


X  drag:  1  cos#  |  +Aq  \  sinf  |)?7  1  17  j  cosi.  (2.24) 

Given  $,  then  Equation  2.21  may  be  used  to  find  7i,  and  the  static  cable  equations  can  be 
integrated  up  to  the  water  surface.  This  gives  the  total  static  configuration. 

The  towfish  dynamics,  imder  the  assumption  of  small  dyneimic  angles,  is  given  by  the 
frequency- domain  equations 


{—{Mp  -1-  Afop)u>®  ut>Bp)P 

=  fj(4  -  ^i)  -  Wcos§4, 

(2.25) 

{—{Mq  -I-  ilfaQ)u>^  -f  iuBQ)Q 

=  fi -Wsinii, 

(2.26) 

(- Jw*  -f  iuBt  -1-  fiAe)i 

=  TjA^i. 

(2.27) 
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where,  as  for  the  cable,  the  damping  is  based  on  describing  functions,  and  to  be  found 
iteratively.  Thus,  we  can  easily  write 


(2.28) 


where  H  is  a  3  x  3  matrix.  We  also  need  a  kinematic  relationship  between  the  motions  of 
the  towhsh  2uid  the  motions  of  the  connection  point: 


p 

Q 

=  G 

fi 

k 

Pi 

9i 


0  1  0 
-10^ 


P 

Q 


:=R 


P 

Q 

$ 


(2.29) 


P 

• 

Q 

Pn 

9n 

Breaking  the  matrix  /3  of  Equation  2.20  into  2x2  blocks  numbered  ySn,  y92i>  ^d 
the  relationship  between  the  towfish  motions  and  the  motions  of  the  top  point  is 


(2.30) 


This  calculation  complete,  one  may  then  use  Equations  2.28  and  2.29  to  get  yj.  The  local 
“state  vectors”  yi  may  then  be  found  all  the  way  through  the  cable  using  Equation  2.18. 

In  order  to  carry  out  iteration  for  the  drag  nonlinearity,  one  may  assume  ^ln  initial 
spatially- varying  damping  &q(s),  6p(s)  as  well  as  initial  values  for  Bp  and  Bq,  and  find  all 
the  vectors  yi.  The  equivalent  linearization  is  then  found  at  each  point;  this  leads  to  a 
different  set  of  damping  coefficients,  which  in  the  next  iteration  gives  a  new  set  of  motions. 
The  iterations  are  continued  until  the  motions  settle  sufficiently. 

This  straight  iterative  loop  has  been  foimd  to  be  stable  in  general,  although  it  may  take 
a  very  long  time  to  converge.  To  speed  up  the  calculations,  a  Newton-Raphson  or  bisection 
approach  may  be  helpful.  For  the  present  work,  bisection  of  the  successive  motion  estimates 
has  been  used,  and  acceptable  results  can  be  obtained  in  less  than  five  iterations,  for  most 
cases. 
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2.4.3  The  Effects  of  Current 


Let  v{s,t)  =  ({s)3in{\t),  and  the  current  be  U.  In  actuality,  v{s,t)  has  some  phase  angle 
associated  with  it,  and  U  may  be  a  function  of  depth,  but  we  forego  these  dependencies 
(without  loss  of  gener^dity)  to  keep  the  following  formulations  simple.  We  investigate  three 
cases  of  quadratic  drag: 

1.  U  =  0‘,  see  Figure  2-5(a).  We  can  use  ■^sin(Xt)  as  the  describing  function  for 
stn(At)  I  stn(At)  |,  so  that 


t»(s,  t)  I  «(s,  t)  |~  C(«)  I  <i»)  I  (2.31) 

The  Foiirier  component  associated  with  this  force  is  therefore  b\{3)  :=  ^C(^)  I  C(^)  I- 

2.  17  >  0,  without  loss  of  generality,  and  U  is  large  enough  that  U  +  C(«)  >  0  everywhere 
on  the  length.  See  Figure  2-5(b).  Then  (t7  +  u)  ]  17  +  r  |  is  positive  everywhere  also. 
This  signal  can  be  decomposed  into  its  Fourier  components;  the  D.C.  component  is 
positive,  and  gives  rise  to  an  offset  in  the  static  shape.  The  first  harmonic  component 
of  (u  +  17)  I  u  +  17  I  is  given  by 

1 

^i(^)  ~  f  J  +  C(^)«t”(.^0  1  V  +  C(s)sin{At)  I  sm(At)dt,  (2.32) 

where  2T  is  the  period  of  oscillation.  Because  17  +  r  is  positive,  we  may  as  well 
integrate  with  the  quadratic  part  replaced  by  a  square.  Applying  some  standard 
integral  identities,  it  is  foimd  that  6f(s)  :=  2UC{s).  This  leads  to  the  approximation 

(t;(s,t)  +  17)  I  v{3,t)  +  17  j~  217C(s)sm(At).  (2.33) 

This  term  reflects  that  fact  that  (r  +  17)(u  +  17)  =  v*  +  217v  +  17®  is  dominated  by 
217 V  when  U  is  large.  The  17®  term  is  canceled  by  the  static  solution. 

3.  17  >  0,  but  there  are  regions  in  each  cycle  where  17  +  w(s,t)  goes  negative.  See 
Figure  2-5(c).  The  crossover  points  /3,(^)  (there  are  two)  can  be  found  in  the  time 
range  {-T,T)  by  setting 
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a)  Case  1:  U  =  0 


time 


time 


Figure  2-5:  Describing  functions  for  three  drag  regimes  based  on  equalized  first  Fourier 
coefficients. 
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C(*)«n  +  tr  =  0  ,  or 


A(-)  = 


(2.34) 


Through  careful  manipulation  of  integrals,  one  is  able  to  deduce  the  first  harmonic 
amplitude  [97]; 


6?(s)  :=  jffn(((j))- 


(2.35) 


where  sffn(-)  is  the  sign  of  the  argument.  It  may  be  verified  that  in  the  limit  as 
C(s)  -*  U,  bf  approaches  bj.  In  addition,  as  U  0,  b^  approaches  b\. 


In  Figure  2*6  are  plotted  the  contours  of  b\  for  a  range  of  values  of  C(3)  aad  of  U .  The 
portion  above  the  U  =  C(^)  line  is  due  to  Case  2  above,  the  portion  on  the  C(«)-axis  is  from 
Case  1,  and  the  area  connecting  the  U  =  C(s)  line  to  the  ^(s)-axis  is  from  Equation  2.35. 
It  can  be  seen  that  bi  depends  quadratically  upon  ^(s)  when  there  is  no  ciurrent,  but  that 
when  the  current  is  large,  this  relationship  is  linear.  Case  3  considerations  provide  a  smooth 
contour  between  these  two  simpler  regimes. 


2.4.4  Harmonic  Balance  Results 

The  data  of  Figtire  2-3  has  been  plotted  against  the  predictions  of  the  harmonic  balance 
method,  in  Figure  2-7.  Each  cmve  here  represents  one  frequency  response  curve  for  a 
particular  amplitude,  each  of  which  is  listed  to  the  left  of  the  curves.  Overall,  the  matching 
between  simulation  results  2uid  the  harmonic  balance  method  is  extremely  good,  although  in 
the  high-amplitude,  high-frequency  regime,  one  notes  some  discrepancy.  This  is  due  to  an  a 
priori  tolerance  specification  in  the  harmonic  balance  algorithm.  An  important  observation 
is  that  the  simulation  required  fourteen  hours  to  run  while  the  harmonic  balance  data  were 
generated  in  approximately  fifteen  minutes,  all  on  a  one  megaflop  workstation. 
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Figtire  2-7:  Comparison  of  simulation  and  harmonic  balance  frequency  response  estimates 
gain  (top)  and  phase  (bottom). 


In  order  to  assess  the  impact  of  ctirrent  on  the  lateral  frequency  response,  the  above 
calc\ilations  were  repeated,  with  a  current  of  two  centimeters  per  second.  As  might  be  antic¬ 
ipated,  the  current  introduces  significant  attenuation  and  phase  loss  in  the  plant;  Figure  2-8 
shows  this  trend  clearly,  when  compared  to  Figure  2-7.  We  see  that  the  curves  associated 
with  the  four  lowest  excitation  amplitudes  are  coincident  because  of  the  linear  drag  of  Case 
2.  The  remainder  of  the  curves  are  in  the  regime  of  Case  3,  which  dictates  that  the  cur¬ 
rent  has  a  diminishing  effect  with  growing  amplitudes.  In  general,  as  the  current  becomes 
larger,  one  sees  more  and  more  of  the  response  controlled  by  linear  damping,  and  a  smooth 
transition  into  amplitudes  that  are  large  enough  that  quadratic  drag  is  still  dominant. 

Since  the  harmonic  balance  method  incorporates  axial  dynamics  with  geometric  coupling 
and  elongation,  calculations  were  also  made  for  the  nominal  system’s  heave  response  and 
coupled  dynamics.  For  the  pure  heave  frequency  response,  wave  periods  between  two  and 
twenty  seconds  were  considered,  with  towpoint  motion  amplitudes  ranging  from  one  to  four 
meters.  The  results  are  shown  in  Figures  2-9  and  2-10;  the  first  shows  the  dependence 
of  towfish  heave  amplitude  on  frequencies,  while  the  second  indicates  the  dynamic  tension 
amplitudes  at  the  top  and  bottom  ends  of  the  cable.  For  this  length  of  cable,  there  are  no 
pronounced  resonance  effects,  as  is  clear  from  the  response  amplitudes.  However,  when  the 
excitation  amplitude  is  large,  the  dynamic  tensions  at  both  cable  ends  can  be  large  also. 
The  two-meter  excitation  curve,  for  exzunple,  suggests  that  snap  loading  at  the  towfish  may 
occur  for  frequencies  above  .8  radians  per  second,  2uid  snap  loading  at  the  surface  may  occur 
above  1.5  radians  per  second.  The  designer  would  be  well-advised  to  modify  the  system 
parameters  in  this  situation. 

The  harmonic  balance  method  was  developed  primarily  for  studying  the  dynamics  of 
mooring  lines,  in  which  significant  curvature  and  geometric  coupling  may  take  place  due  to 
gravity  alone.  For  a  heavy,  towed  vehicle  system,  we  have  seen  that  the  lateral  response 
of  the  cable  in  the  wave  frequency  range  is  confined  to  a  segment  near  the  surface.  Thus, 
when  both  surge  and  heave  motions  are  imposed  at  the  surface,  it  is  generally  the  heave 
alone  which  is  observed  at  the  towfish.  This  point  is  illustrated  in  in  Figure  2-11,  where 
the  vertical  and  horizontal  towfish  responses  are  given,  for  a  towpoint  motion  of  two  meters 
surge  and  heave,  with  the  heave  leading  the  surge  by  ninety  degrees.  The  current  for  this 
run  is  one  meter  per  second,  which  leads  to  a  static  towfish  offset  of  229  meters.  The 
dynamic  tension  curves  are  indistinguishable  from  those  in  Figure  2-10,  which,  along  with 
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Figure  2-8:  Harmonic  balance  frequency  response  estimates  with  a  current  of  .02  me- 
teis/second:  gain  (top)  and  phase  (bottom). 
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Figure  2-9:  Amplitude  of  towfish  heave  frequency  response. 


frequency,  rad/sec 


Figure  2-10:  Cable  tension  frequency  response  during  heaving. 
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Figure  2-11:  Surge  and  heave  towiish  response  during  coupled  motions. 


Figure  2-11,  leads  one  to  conclude  that  the  effects  of  geometric  coupling  are  small  for  heavy 
cable  systems  towed  at  low  speeds.  One  must  remember,  however,  that  the  static  tension 
induced  by  towing  may  be  significantly  larger  than  that  due  to  weight  alone;  hence,  the 
total  tension  can  be  affected  by  towing. 


2.5  Perturbation  Analysis 

Our  second  approximation  of  the  lateral  dynamics  of  a  towed  cable  is  a  more  classical  per¬ 
turbation  study  using  the  multiple-scale  technique,  with  the  first-order  in-plane  equations 
(to  be  detailed).  As  perturbation  methods  work  best  near  natural  modes,  the  generality  of 
harmonic  balances  is  lost.  Nonetheless,  the  perturbation  approach  is  worthwhile  because  it 
yields  simple  algebraic  relationships  that  govern  the  response,  and  provides  a  closed-form 
solution  for  limit  cycling  behavior  in  the  closed  loop. 

2.5.1  First-Order  Three-Dimensional  Equations  for  a  Vertical  Cable 

When  the  cable  is  known  to  be  nearly  vertical,  an  expansion  may  be  used  which  is  more 
suitable  than  that  for  the  equivalent  linearization  approach  above.  In  particular,  we  now 
allow  for  higher-order  terms,  to  be  indexed  by  the  small  parameter  c  which  characterizes 
the  excitation  amplitude.  It  is  known  from  physics  that  the  tension  T  and  axial  velocity 
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u  contain  only  even-powered  terms  in  c,  and  that  v,  to,  d  all  contain  only  odd-powered 
terms  in  Thus, 

T  = 

<t>  =  €<j>\  -i-  -I-  0(€*), 

6  =  €$1  -f  €*^3  0(£®), 

«  =  tio  -f-  €*tt2  +  0(«^), 
t>  =  €Vi  -I-  «*t>3  +  0(c®),  and 
to  =  €tOi  -I- €®t03 -I- 0(e®). 

Although  we  have  dropped  the  arguments,  it  must  be  remembered  that  the  tension,  angles, 
and  velocities  are  functions  of  t  and  s.  In  addition,  since  we  will  not  be  discretizing  the  cable, 
the  subscripts  of  variables  in  this  section  refer  to  the  associated  order  in  the  expansion,  not 
the  spatial  coordinate  as  before.  For  our  physical  scenario,  the  static  solution  (denoted  by 
an  “0”  subscript)  is  given  by: 

To  =  W  +  w„3,  (2.36) 

uo  =  0, 

where  the  coordinate  s  is  assiuned  zero  at  the  bottom  of  the  cable.  We  see  that  the  zeroth- 
order  solution  is  one  in  which  there  is  no  motion,  and  the  tension  is  created  only  by  the 
weight  of  the  system. 

We  now  write  out  the  three-dimensional  equations  in  the  first-order  approximation:  they 
are 


=  EI^-To^  +  7:R3{4>u0uVi,w^),  (2.37) 

'If  the  top  of  the  cable  is  excited  with  a  lateral  in-plane  motion  (sm(At),  for  example,  T  and  u  respond 
the  same  irrespective  of  the  sign  of  c.  On  the  other  hand,  the  lateral  motions  clearly  depend  on  e’s  sign. 
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dvi  _ 
ds  ~  at  ’ 

dwi  _  601 

ds  ~  dt' 

It  is  important  to  note  that,  except  for  drag  forces,  these  equations  are  uncoupled;  that  is, 
two  of  the  equations  describe  purely  in-plane  motions,  while  the  other  two  describe  purely 
out-of-plane  motions.  If  q  and  r  are  in-plane  and  out-of-plane  deflections,  respectively,  then 
the  relevant  equations  are: 


(m-f  maj— 
(m  +  rria)-^ 


=  T 


d^q 

ds^ 


+  (2.J8) 


~  dt^ [dtj  ■‘■(at) 


In  the  presence  of  ambient  current,  the  normal  component  should  be  added  to  ^  above, 
and  the  binormal  component  to 

We  do  not  consider  higher-order  expansions  of  Equations  2.12,  as  they  have  been  carried 
out  by  Howell  [56]  to  investigate  the  geometric  nonlinearities.  More  importantly,  we  do  not 
anticipate  that  the  coupling  which  comes  with  the  higher-order  approximations  is  necessary 
for  our  scenario.  The  motions  are  assumed  to  be  very  slow  and  overdamped,  and  the  cable 
system  is  quite  heavy,  maintuning  its  nearly  vertical  orientation. 


2.5.2  First-Order  Two-Dimensional  Equations  for  a  Vertical  Cable 

When  only  the  in-plane  motions  are  considered,  the  zeroth-order  solution  is  the  same  as  in 
the  three-dimensional  case,  and  only  two  first-order  equations  need  to  be  considered: 


/  ,  .dvi 

(m  -I- 


dt 


^  d^6i 

—  El  -I-  SH2(^i>0,vi,0), 


dVi 

17’ 


Equation  2.38  can  be  used  for  the  evolution  of  the  lateral  displacement  q. 


(2.39) 
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2.5.3  Linear  Analytical  Solution  for  In-Plane  Motions 

There  exists  an  exact  analytical  solution  to  the  first-order  problem  in  the  two-dimensional 
case  when  there  is  no  bending  stiffiiess  and  the  damping  is  linear.  This  linear  solution  is 
well  known;  however,  working  it  out  in  this  section  helps  to  form  a  foundation  for  our  later 
results. 

From  Equations  2.40  and  2.21,  the  first-order  equations  of  motion  for  the  cable,  the 
towfish,  and  the  compatibility  relation  are  given  respectively  by: 

{w„s  +  -t-  -  bvi, 

(2.40) 

dt  ’ 


,  ,  y,dVl 

dv\ 

ds 


where  b  and  B  are  the  linear  damping  coefficients  of  the  cable  and  towfish,  respectively. 
Note  that  an  tmderlined  quantity  (e.g.,  Vi)  indicates  it  is  to  be  evaluated  at  the  lower  end 
of  the  cable.  Similarly,  a  quantity  with  an  overbar  (e.g.,  vi)  indicates  that  it  is  to  be  taken 
at  the  top  of  the  cable.  We  make  use  of  a  number  of  nondimensional  variables: 


s 


$  :=  <f)—,  and 
d 


Then  the  Equations  2.40  can  be  written  in  nondimensional  terms  as  follows: 


=  (fcio- 

^  ,  and  (2.41) 

dui  _  dOi 
da  dr  ’ 
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where 


*1 

*2 

kz 

k4 

is 


W 

(m  +  m«)s’ 
W 

X(m  +  ma)ff’ 


{Mp  +  Map)g' 


B  fl 

Mp  +  MaP  Y  9 


We  defiiie  a  new  spatial  variable  V*  •=  k\<T  +  ibjt  so  that  taking  the  derivative  of  the  cable 
equation  of  motion  with  respect  to  a  and  using  the  compatibility  relation,  the  following 
result  is  obtained: 


(2.42) 


We  now  use  separation  of  variables  with 


^i(V’>t)  =  R{xl’)NiT)  =  a(e’^’’  +  cc)i?(^), 


(2.43) 


where  cc  indicates  the  complex  conjugate  of  the  preceding  term,  a  is  a  scalar  amplitude, 
and  A  is  a  nondimensional  excitation  frequency.  With  the  definition 


a  := 


A2  -  ikzX 


a  form  of  Bessel’s  equation  is  obtained  for  the  first  phasor  component: 


,<PR 


dR 


V*  +  2V'3-r  +  i/^aR  =  0 


(2.44) 


Since  a  is  complex,  it  is  helpful  to  make  a  change  of  variables  before  giving  the  solution. 
Letting  the  complex  spatial  variable  z  be 


2  :=  2\/a^, 
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and  performing  some  algebra  gives  the  solution  of  Equation  2.44  as 

W)=  +  (2.45) 

where  Jk  is  the  ibth-order  Bessel  function  of  the  first  kind  and  Yk  is  the  kth-order  Bessel 
function  of  the  second  kind.  The  scalar  c  is  to  be  determined  from  the  lower  boundary 
condition,  and  the  coefficient  of  Ji(2)  has  been  absorbed  into  a.  This  leads  to  an  expression 
for  the  velocity  in  the  cable: 


The  constsmt  c  may  be  found  easily  from  the  bottom  boundary  condition  as 


(2.46) 


c 


+  zaJojl) 

k4Yiiz)+^aY„{z) 


,  where 


(2.47) 


-A*  +  ih\ 
- 20*1 


Finally,  the  amplitude  a  may  be  specified  from  the  top  boimdary  condition.  We  impose 


i/  =  t'Ac’*’’  +  cc. 


(2.48) 


which  leads  to  the  final  unknown  for  this  problem: 

_  — 2a*i 

Jo(z) +  €?;(«)' 

The  amplitude  and  phase  response  of  the  system  may  be  found  from  these  relations,  by 
taking  the  magnitude  and  the  argument  of  a.  When  the  mass  and  damping  characteristics 
of  the  clump  weight  go  to  zero,  one  requires  that  c  =  0,  leading  to  the  classical  shape  for  a 
freely  hanging  cable  [19].  For  the  case  with  nonzero  towfish  mass  and  weight,  an  example 
Bode  plot  is  given  in  Figure  2-12.  The  response  is  typified  by  an  infinite  number  of  peaks, 
and  the  figure  shows  that,  as  for  finite-dimensional  systems,  added  damping  shifts  the  peaks 
to  lower  frequencies  with  lower  gain.  The  first  three  natural  frequencies  lor  the  nominal 
system  are  0.0982,  0.2678,  and  0.4725  radians  per  second. 
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degrees 


Figure  2-12:  Bode  plot  for  hanging  cable  with  three  values  of  linear  damping. 
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The  control  system  designer  must  be  aware  that  in  applications  where  the  plant  has 
minimal  damping,  the  presence  of  multiple  peaks  above  the  primary  cutoff  frequency  may 
impose  restrictive  performance  conditions.  Furthermore,  the  presence  of  currents  in  the 
fluid  has  no  effect  on  the  response  with  the  linear  drag  law. 

2.5.4  Nonlinear  Frequency  Response 

In  this  section,  we  try  to  xmcover  more  rigorously  the  mechanism  that  leads  to  the  amplitude- 
varying  frequency  response  in  Figure  2-3.  The  goal  of  the  analysis  is  to  characterize  the 
shape  of  the  figure  using  a  few  algebraic  relations,  instead  of  through  simulation  or  harmonic 
balance.  An  important  simplification  is  that  we  will  neglect  nonlinear  drag  on  the  towfish 
itself.  This  condition  is  reasonable  given  the  very  large  area  of  the  cable  compared  to  the 
towfish,  and  is  necessary  for  the  analysis  to  proceed  smoothly. 

We  begin  with  the  same  notation  as  above.  Here,  however,  the  response  is  assumed 
to  be  a  perturbation  of  one  of  the  natural  modes.  As  such,  the  tmdamped  natural  modes 
first  must  be  identified  by  setting  the  motions  at  the  upper  boimdary  equal  to  zero.  The 
variables  kz  and  ks,  which  quantify  the  eflects  of  linear  drag,  are  set  to  zero,  and  quadratic 
drag  is  incorporated  in  the  cable  equation  as  follows: 

di/  do 

~  =  {kx<T  -I-  *2)^  kx9  -  €{v -It  u)  \  V  +  u  I,  (2.50) 

where  c  =  2^^^)  •  should  be  noted  that  this  is  a  new  expansion  in  that  vx  and  $x 
of  Equation  2.42  have  been  replaced  by  v  and  9.  In  other  words,  we  have  adopted  what 
was  the  first-order  equation  for  the  physical  cable  to  describe  the  complete  system  for  this 
analysis.  For  applications,  e  generally  will  be  between  .1  and  .3,  so  a  perturbation  solution 
of  Equation  2.50  can  be  expected  to  give  reasonably  good  resxilts.  Also,  the  effects  of  a 
nondimensionalized  ambient  current  u  have  been  included  in  the  cable  equation. 

To  carry  out  multiple-scale  an2Llysis  (e.g.,  [16]),  we  let  r  be  a  slow  time  scale  such  that 
r  =  e^r.  The  expansion  is  only  to  be  csurried  to  third-order  in  e,  so  only  one  slow  time  scale 
is  necessary.  The  variables  v  and  6  thus  become  functions  of  z,  r,  and  r,  and  they  may  be 
expanded  in  terms  of  e  as  follows: 
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V 


(2.51) 

(2.52) 


=  €vi  +  +  J<t)  +  0(€*), 

^  =  c^i  +  €*(^3  +  ^)  +  0{e^), 

where  t)  is  the  nondimensional  cable  displacement.  The  top  excitation  should  be  O(e^);  it 
can  be  written  as: 


fj  =  +  cc).  (2.53) 

The  scalar  h  dictates  the  nondimensional  size  of  this  excitation. 

The  third-order  equation  comes  out  as  follows: 


—  {kic  +  ^2)-^  —  ^1^3  =  — ^  +  h{a>^  +  Ai)(e’*'^  -f-  cc)  -  (i/j  -|-  u)  |  i/j  u  |  .  (2.54) 

Secular  terms  indeed  appear  in  the  right-hand  side,  so  we  can  pursue  the  perturbation 
solution  by  setting  them  equal  to  zero. 

Toward  this  end,  the  first  two  terms  on  the  right-hand  side  of  Equation  2.54  are  straight¬ 
forward,  but  the  l2ist  must  be  treated  with  the  same  care  as  in  the  harmonic  balance.  We 
use  precisely  the  same  results  here  as  given  in  Section  2.4.3,  with  the  following  link.  Suppose 
first  that  the  amplitude  a  is  replaced  by  a  slowly- time- varying  equivalent: 


Mr)  =  G(T)e^(^, 


where  Cr(T)  ^lnd  f (t)  are  real.  Then  defining 


(2.55) 


<{r,  2) 


<^nGiT) 

aki 


where  J{z)  :=  Jo{z)  -i-  cYo{z)  and  u>„  is  a  natural  frequency  of  the  system,  leads  to 


»'i(t,t,  2)  =  C(r»^)W‘*'n’‘  +  ^(r)).  (2.56) 

The  complete  right-hand  side  of  Equation  2.54,  for  the  three  drag  cases  of  Section  2.4.3, 
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is  as  follows: 


Case  1:  |  (  |  +  cc) 

^^(xe*"’*’’^+cc)7(2)+/i(<rA*+*i)(ie‘^’‘+cc)+<  Case  2:  t/^C(te^«'^e^(^)  +  cc) 

Case  3:  +  cc)  J 

(2.57) 

Recalling  that  <7  =  we  make  the  following  definitions  so  that  the  Galerkin  proce¬ 

dure  can  be  used  to  eliminate  the  spatial  dependence: 


Qi  := 


03  := 


05(1)  := 


=  j  J^{z)dz 
=  J  z^J{z)dz 
=  j^J{z)dz 

~  Iz  I  I 

=  h\{T^z)J{z)dz^ 


76(r)  := 


h\^a2  *2tx2 

kihaz 

4w^Q4 

Uf^nai 

aki 

Qs(z) 

2  ■ 


Making  the  substitutions  into  Equation  2.57  and  dividing  through  by  the  following 

requirement  is  obtained: 
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IaG  I  G  1 1 

HsGi 

76* 


^=0,  (2.58) 


7i*(^  +  Gi^)  +  (72  +  73)(co6(tfr  -  ^)  +  isin(ST  -  f))  + 


where  6  :=  ^  Following  Nayfeh  and  Mook  [74]  in  expressing  the  steady-state  amplitude 
as  a  function  of  the  detiming  parameter  we  look  at  the  three  drag  cases  separately: 


1. 


-7i^*  +  +  474(72  +  73)’  /72  +  73 


27I 


74 


when  ^  »  0.  (2.59) 


2. 


+  7I 


(2.60) 


3.  The  equation 


7l  +  7?^'<?'-(72  +  73)'  =  0  (2.61) 

must  be  solved.  Since  the  amplitude  of  this  case  is  known  to  be  less  than  that  of  Case 
1,  it  suffices  to  decrement  from  the  Case  1  result  imtil  Equation  2.61  is  satisfied. 

2.5.5  Trends  of  the  Nonlinear  Frequency  Response 

An  investigation  of  the  above  equations  indicates  that  there  are  several  simple  mechanisms 
at  work.  First,  the  analysis  always  predicts  attenuation  away  from  reson2mt  frequencies 
through  Equations  2.59  and  2.60.  This  trend  is  coincident  with  intuition  and  with  simu¬ 
lation  results,  as  shown  in  Figure  2-7.  However,  the  perturbation  technique  is  only  valid 
near  a  natural  frequency;  Figure  2-13  shows  how  the  analysis  compares  with  simulation  re¬ 
sults,  in  complete  analogue  with  Figure  2-7.  The  perturbation  solution  was  computed  with 
respect  to  the  three  natural  modes  in  the  frequency  range,  switching  at  the  discontinuities 
in  the  figure.  One  finds  that  the  matching  is  very  good  in  the  immediate  neighborhood  of 
the  natur£d  modes,  but  is  quite  poor  at  any  distance  from  them.  Presumably,  the  fact  that 
the  aescribing  functions  are  always  applied  to  the  natural  mode  shapes,  even  away  from  the 
natural  frequencies,  is  responsible  for  the  poor  fit. 
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Figure  2-13:  Comparison  of  nonlinear  frequency  response:  simulation  vs.  perturbation 
approximation. 


An  important  second  characteristic  of  this  perturbation  analysis  is  that  |  (?  |a  so 
the  gain  is  inversely  related  to  the  square  root  of  the  amplitude.  Equation  2.59  indicates 
this  effect,  which  is  roughly  consistent  with  the  simulation  results  of  Figure  2-13,  especially 
near  the  natural  modes.  Thus,  in  one  equation,  we  have  characterized  the  mechanism  that 
makes  the  filteiing  action  of  the  cable  stronger  for  large  velocities. 

With  respect  to  current  effects,  it  is  clear  from  Figure  2-6  that  the  6i-contours  are  forced 
leftward  in  significant  currents.  This  fact  shows  explicitly  that  the  frequency  response  is 
reduced  in  the  presence  of  current,  consistent  with  the  harmonic  balance  predictions.  In 
fact,  increasing  the  currents  makes  the  sustainable  response  (  approach  zero  asymptotically, 
as  for  a  given  amplitude  of  oscillation  may  be  arbitrarily  large.  It  is  also  interesting  to 
note  that  the  maximum  sensitivity  of  the  bi  contours  to  changes  in  current  is  along  the 
U  =  C  line;  the  sensitivity  is  very  low  for  small  currents. 

Figure  2-14  shows  a  comparison  of  the  frequency  response  in  current  found  by  sim\ilation, 
harmonic  balance,  and  the  perturbation  method.  The  excitation  for  which  these  calculations 
were  made  is  denoted  by  the  “x”  in  Figure  2-13  (.093  radians  per  second,  2.154  meters 
amplitude).  The  performance  of  the  harmonic  balance  technique  is  good,  except  that  the 
simulations  indicate  a  “saturation  current,”  beyond  which  there  is  no  measured  response. 
The  towfish  surge  response  predicted  by  harmonic  balance  smoothly  approaches  zero  as  the 
current  is  increased.  The  m\iltiple-scale  prediction  is  accurate  when  the  cxirrent  is  small, 
but  becomes  increasingly  poor  as  the  current  is  increased,  and  again  the  complete  loss  of 
the  response  above  the  saturation  is  not  evident.^  The  deterioration  of  accuracy  in  large 
currents  is  not  unexpected,  given  the  large  amount  of  approximation  involved  and  the  fact 
that  the  cable  motions  may  actually  have  very  little  similarity  to  any  natural  modes. 

2.6  Control  Considerations 

2.6.1  Model  Reduction 

The  models  described  so  far  have  used  a  number  of  assiunptions  about  the  geometry  (e.g., 
first-order  equations  only)  and  the  physics  (form  of  the  drag  law,  etc.).  The  transfer  function 
from  top  motions  to  bottom  motions  is,  however,  still  infinite-dimensional,  as  the  usual 

’The  steps  in  the  anajytical  approximation  are  a  result  of  the  incremental  procedure  used  to  solve  Equa¬ 
tion  2.61. 
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Figure  2-14:  Comparison  of  attenuation  of  frequency  responses  in  current:  simulation, 
harmonic  balance,  and  multiple-scale  analysis. 
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numerical  procedures  for  modeling  this  type  of  system  may  have  fifty  or  more  nodes. 

There  is  available  a  large  body  of  literature  on  the  feedback  control  of  linear  infinite¬ 
dimensional  systems  (e.g.,  [31],  [28]).  The  main  issues  are  quantification  of  the  effects  of 
finite-dimensional  control,  and  alternatively,  the  reduction  of  plant  model  order  while  main¬ 
taining  small  modeling  errors  [69],  [85].  These  techniques  are  largely  based  on  frequency- 
domain  properties,  and  there  are  by  comparison  very  few  tools  av^able  for  general  nonlinear 
problems. 

For  this  thesis,  model  reduction  is  carried  out  using  an  ad  hoc  time-domain  technique 
[55].  Essentially,  the  reduction  process  is  viewed  as  a  system  identification  in  which  the 
structure  of  the  resulting  model  is  highly  constrained;  in  fact,  it  is  a  low-order  finite- 
difference  approximation  (see  Section  2.7).  Data  from  either  physical  tests,  or  from  high- 
order  models  are  taken  as  “truth,”  and  an  iterative  se2irch  in  the  parameter  space  is  con¬ 
ducted  to  minimize  the  simulation  errors  of  the  reduced  model  [18].  As  this  approach  has  no 
guarantees,  the  best  we  can  do  is  carry  out  the  identification  for  a  truth  I/O  data  set  that 
has  the  expected  characteristics  of  the  closed-loop  maneuvers  embedded  in  it.  Modeling 
error  issues  are  revisited  in  Chapters  4  and  5. 

The  usual  manifestions  of  model  reduction  errors  in  feedback  control  are  observation 
and  control  “spillover”  [6].  Control  spillover  refers  to  the  excitation  of  high-frequency  plant 
modes  by  the  controller,  and  observation  spillover  refers  to  the  contamination  of  state 
measurements  by  these  modes.  Observation  spillover  can  be  reduced  greatly  by  appropriate 
filtering,  so  closed-loop  stability  is  not  usually  a  problem  when  the  plant  is  stable.  The 
performance  of  such  a  feedback  system  may  suffer,  however,  due  to  control  spillover,  and 
the  situation  is  exacerbated  when  the  plant  is  lightly  damped.  The  next  two  sections  suggest 
that  the  cable  of  interest  will  not  siiffer  large,  high-frequency  motions. 

2.6.2  Energy  Dissipation  in  Vibratory  Modes 

For  a  long  cable  in  water,  it  can  be  argued  that  the  effects  of  high-frequency  excitation 
are  negligible,  due  to  the  very  large  drag  terms  associated  with  such  motions.  We  consider 
a  cable  under  constant  tension  with  fixed  ends,  in  a  viscous  fluid.  The  trivial  eigenvalue 
problem  gives  the  n-th  natural  frequency  as  follows: 


nr 

1  ^ 

1  m  +  rria 
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We  wish  to  characterize  how  energy  is  dissipated  into  the  fluid,  and  one  useful  measure  is 
the  ratio  of  the  energy  dissipated  in  a  quarter-cycle  to  the  initial  kinetic  energy  associated 
with  the  mode.  With  q  as  the  amplitude  of  the  mode  shape,  the  kinetic  energy  turns  out 
to  be 


AL  ' 


(2.62) 


To  get  the  energy  lost  to  the  fluid  in  one  quarter  cycle,  we  must  evaluate  the  integral 
El  =  sin*(^^)coj*(a»„t)  |  stn(^^)cos(w„t)  |  dadt, 

where  v  =  u>n9.  Because  the  modes  are  symmetric  in  space  and  because  cos(u7„t)  >  0,  t  G 
[0,Tr/2],  we  have 


El  =  -pCddv^n  /'  sin^(^^~)cos^{unt)d3dt. 

2  Jo  Jq  L 

Separating  the  integrals  and  using  standard  identities,  we  And 


(2.63) 


„  ApCddq^Tri^ic 

^^  =  '9(rTrr^)i’” 

El  _  GApCjq 
Ek  ~  9v^d{p^  Pc)' 

For  a  typical  deep-water  system,  this  ratio  is  on  the  order  of  lOg,  for  steel  cables,  and  50g  for 
synthetic  cables  (with  g  expressed  in  meters).  Inspection  shows  that  a  larger  diameter  and 
cable  density  leads  to  a  smaller  ratio;  that  is,  there  will  likely  be  more  residual  oscillations. 
For  modal  amplitudes  greater  than  a  few  centimeters,  however,  the  initial  kinetic  energy  is 
generally  lost  well  before  the  quarter-cycle  is  completed,  suggesting  that  no  overshoot  will 
occm  in  the  real  system. 

Very  small-scale  motions  that  do  allow  oscillation  may  correspond  to  a  very  low  Reynolds 
number,  for  which  the  asstimed  quadratic  law  law  should  be  replaced  by  a  linear  one.  This 
change  can  be  incorporated  for  Re  <  1000  [81].  Calculations  show  that  when  oscillatory 
motions  are  admitted  by  the  energy  argument,  the  linear  damping  ratio  comes  out  to  be 
approximately  0.9. 
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2.6.3  Limit  Cycling  in  the  Closed  Loop 

Another  way  to  quantify  the  effects  of  modd  reduction  is  to  analyse  the  limit'cycling  behav¬ 
ior  in  closed-loop  operation.  It  often  happens  in  applications  that  the  robustness  conditions 
cannot  be  met  at  all  points  in  the  state  space.  For  example,  the  natural  damping  in  a 
towed  cable  may  tend  to  sero  as  the  operating  speed  is  reduced,  and  the  standard  margins 
supplied  with  an  LQ  regulator  might  well  be  exceeded.  In  such  a  case,  however,  the  system 
does  not  go  wildly  unstable;  what  we  usually  see  is  the  state  entering  into  a  limit  cycle. 
Intuitively,  the  instability  mechanism  is  hdd  in  check  by  the  growing  drag  on  the  cable  as 
the  motions  get  larger. 

The  describing  function  technique  [45]  has  historically  been  the  most  common  way  to 
estimate  analytically  the  size  of  a  limit  cycle.  In  this  procedure,  only  linearized  plants  are 
considered,  and  one  looks  for  the  linearizations  which  cause  the  Nyquist  plot  to  cross  the 
critical  point  (-1).  These  linearizations  then  serve  as  estimates  of  the  states  associated 
with  the  stability  boundary.  The  method  is  an  approximate  one,  and  offers  no  guarantees; 
one  need  only  look  at  the  many  counterexamples  to  Aizermaim’s  conjecture  (see,  e.g.  [92], 
11021). 

The  alternative  route  we  pursue  in  this  section  is  a  perturbation  study  of  the  closed-loop 
system.  The  approach  is  inherently  the  same  as  the  describing  function  technique  in  how 
it  linearizes  the  system,  and  thus  is  subject  to  the  same  pitfalls.  However,  by  extending 
our  earlier  perturbation  technique,  we  can  easily  consider  the  full  infinite-dimensional  plant 
under  linear  feedback  control  and  arrive  at  simple  algebradc  equations  that  govern  the  closed- 
loop  response.  The  describing  function  technique,  on  the  other  hand,  requires  evaluation 
of  a  large  number  of  linearized  plants.  As  before,  we  expect  our  approximation  to  be  good 
for  frequencies  near  the  natural  modes  of  the  plant,  with  small  e. 

We  build  upon  the  conventions  of  Section  2.5.4,listing  only  the  terms  that  come  out 
differently  for  the  current  situation.  The  control  law  is  of  the  form: 

q  =  -kji,  (2.66) 

implying  through  Equation  2.52  that  the  compensator  transfer  function  is  O(e^). 

With  k  consisting  of  a  real  and  an  imaginary  part,  we  set  7  Redefining 
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7j  := 
7s  := 


we  obtain  a  result  similar  to  that  in  Section  2.5.4: 


7i*(-^  +  +  (72  +  { 


74G  I  G  I  i 

7BGt  =  0.  (2.67) 

76f 

Separating  the  real  and  imaginary  parts,  and  looking  at  the  steady-state  solution  only,  one 
finds  that  the  real  part  of  Equation  2.67  does  not  depend  on  the  drag  term.  We  thus  have 
the  frequency  perturbation 


<mz)  _  IflT  +  73r 
dr  7i  ’ 


(2.68) 


for  aU  three  cases,  where  the  *‘r”  subscript  denotes  the  real  part,  and  **1”  denotes  the 
imaginary  part. 

To  find  the  amplitudes,  the  three  drag  cases  are  as  follows. 

1.  The  imaginary  part  of  Equation  2.67  gives 


7i-^  +  (72*  +  l*i)G  -I-  "jiG  I  G  (=  0.  (2.69) 

From  the  definitions,  it  is  easy  to  see  that  71  >  0,  and  73  >  0,  so  that  the  quadratic 
drag  is  stabilizing.  The  impact  of  the  control  transfer  function  ib  is  in  72^  -I-  74,, 
which  can  easily  be  negative  for  a  range  of  values  of  k,  depending  on  the  mode.  If 
72i  +  74*  >  0,  then  we  predict  no  limit  cycle  for  the  associated  mode;  if  the  converse 
is  true,  then  the  steady-state  amplitude  can  be  written  simply  as: 


I  G  1=  -2;2i±2li.  (2.70) 

73 

2.  In  the  presence  of  a  large  current  which  exceeds  the  local  cable  velocity  at  all  points 
on  the  cable.  Equation  2.68  is  still  valid,  but  the  imaginary  part  of  the  equation  now 
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Tkus,  depending  on  the  control  fc,  we  have  dtha  first-order  instability  or  stability. 
Instability  leads  to  Case  3;  alternatively,  if  Equation  2.71  is  stable,  the  decaying 
amplitude  G  stays  in  the  regime  of  Case  2. 

3.  When  vi  +  u  has  zero  crossings,  matters  are  complicated  significantly,  as  indicated  in 
Section  2.4.3.  The  condition  associated  with  the  imaginary  part  is 

dG 

7i^  +  (72t  +  1ax)G  +  7e  =  0,  (2.72) 

which  has  no  simple  dependence  on  C?,  unlike  Cases  1  and  2.  Here  it  is  necessary  to 
iterate  a  few  times  to  get  an  appropriate  value  for  G  in  the  steady-state. 

The  frequency  deviation  is  expressed  directly  in  terms  of  the  real  part  of  the  compen¬ 
sator,  while  the  amplitude  is  related  to  the  imaginary  part  of  the  compensator.  This  is 
clearly  a  restriction,  since  a  purely  real  compensator  (proportional  feedback)  leads  to  limit¬ 
cycling  behavior.  However,  the  accuracy  of  the  perturbation  result  is  unclear  when  the 
compensator  is  complex,  or  even  purely  imaginary. 

For  the  purposes  of  validation,  we  consider  a  compensator  which  is  of  the  proportional- 
plus-integral  (PI)  type.  This  form  was  chosen  because  it  is  simple,  proper,  and  easily 
programmed,  and  because  the  loop  transfer  function  thus  generated  is  of  the  general  shape 
that  we  would  like  to  have  (see  Figure  2-12).  With  ^  =  kp  +  ^,  the  comer  frequency 
is  given  the  values  0.01,  0.02,  0.05,  0.1  (first  mode  of  the  nominal  system),  0.2,  and  0.5 
radians/second.  By  varying  and  looking  at  these  six  cases,  we  can  investigate  the  effects 
of  varying  control  bandwidths  and  phase  angles  in  the  compensator. 

Figure  2-15  shows  the  results  of  this  analysis,  with  each  simulation  point  drawn  as  an  “x,” 
and  the  approximation  taken  about  the  first  mode  only.  Higher  modes  were  computed,  but 
only  the  odd  modes  had  nonzero  amplitudes,  and  those  were  all  several  orders  of  magnitude 
less  than  that  associated  with  the  first  mode.  For  the  nominal  system,  e  =  .205,  and  since 
I  ib  I  is  assumed  to  be  of  order  e^,  we  would  not  expect  good  results  much  above  kp  =  .04. 


The  figure  indicates  that  this  is  not  the  case,  as  we  obtain  good  accuracy  up  to  kp  =  .1  for 
all  six  of  the  cases  considered. 

The  response  rolloff  predicted  by  the  analysis  can  be  explained  as  follows:  large  values 
of  kp  cause  the  frequency  deviation  to  grow,  through  Equation  2.68  ,  and  thus,  in  the 
control  law  is  reduced.  Equation  2.70  then  predicts  that  the  amplitude  will  be  affected;  in 
the  case  of  the  first  mode,  it  is  reduced. 

This  rolloff  is  in  sharp  contrast  to  the  simulation  residts,  in  which  the  amplitude  of 
the  limit  cycle  is  quite  large  at  ibp  =  1.  fri  fact,  the  system  is  unstable  for  kp  >  .9  when 
the  comer  frequency  is  .2  radians/second,  and  for  kp>  A  when  the  comer  frequency  is  .5 
radians/second.  These  situations  correspond  to  a  "sub-first”  mode,  in  which  the  cable  is 
being  drawn  through  the  water  so  rapidly  that  it  loses  180  degrees  of  phase  below  the  first 
natural  mode.  This  loss  of  phase  is  easily  seen  in  Figures  2-3  and  2-7.  The  sub-iiTSt  mode 
instability  is  an  important  motivation  for  the  nonlinear  controllers  we  pursue  in  Chapter  4. 

Another  observation  is  that  the  amplitudes  reported  by  both  the  analysis  and  simula¬ 
tions  are  quite  small  when  a  reasonable  controller  design  is  used.  This  finding  is  consistent 
with  our  earlier  arguments,  and  suggests  that  well-chosen  feedback  laws  may  be  unstable 
when  applied  to  an  undamped  cable  system,  yet  can  be  quite  acceptable  in  practice.  Fig¬ 
ure  2-15(d),  for  example,  implies  that  even  placing  the  comer  frequency  at  the  first  resonant 
mode,  with  unity  gain  at  that  point,  leads  to  motions  of  under  one-meter  amplitude. 

A  final  note  is  that  attenuation  of  the  limit-cycle  is  predicted  when  current  acts  upon 
the  cable.  For  example,  with  U  =  .01  meters  per  second  and  a  comer  frequency  of  .5  radians 
per  second,  the  analysis  predicts  no  limit  cycle  up  to  Jbp  =  .45,  and  a  reduced  limit  cycle 
above  that  point.  However,  because  this  analysis  does  not  take  into  account  the  loss  of 
plant  phase  which  accompanies  currents,  such  results  are  generally  inaccurate. 

2.7  Lumped-Mass  Models 

In  this  section,  we  describe  otir  standard  lumped-mass  model  for  the  lateral  dynamics  of 
the  cable  system.  This  model  will  be  used  through  the  rest  of  the  thesis.  Beginning  with 
the  drag-coupled  first-order  Equations  2.38,  we  discretize  the  cable  into  segments  of  length 
ds,  and  adopt  the  standard  O(ds^)  centered  finite- difference  approximations  to  the  spatial 
derivatives: 
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9i+l  ~  2ft  +  gt-1 
ds^ 


qi+1  -  g»-i 
2ds 


Making  the  substitutions,  we  have 


-  laiiS  ^  +  'h  (2-73) 


The  top  boundary  condition  of  the  cable  is  an  imposed  displacement  of  the  endpoint.  The 
first  node  in  the  discretization  should  thus  be  placed  a  length  ds  down  from  the  surface,  so 
that  Equation  2.73  is  valid  there.  At  the  lower  boundary  is  a  clump  weight  whose  dynamic 
response  is  given  by 


{M  +  Map\^  =  -  \pCd,tou,fi,hApg^^^  +  rj  (2.74) 

(M  +  MaR)rit  =  Wr,  -  \pCa,u„vfUhARrty/^ 

In  a  model  with  n  nodes,  by  definition  we  impose  the  vehicle  dynamics  on  the  first  node, 
with  the  0(ds)  forward-differences  approximation  =  g,,i  ~  (52  -  9i)/ds. 

This  formulation  does  not  quite  lead  to  a  simple  mass-spring  system;  it  can  be  e2isily 
shown  that  a  carelessly  constructed  mass-spring  approximation  has  a  fund2unentally  differ¬ 
ent  response,  even  £is  the  number  of  nodes  goes  to  infinity.  In  addition,  it  is  noteworthy 
that  linearizations  of  this  model  do  not  have  a  pure  delay,  nor  any  zeros.  The  lack  of  zeroes 
is  consistent  with  Figure  2-12,  and  for  most  systems,  the  time  constant  is  several  orders  of 
magnitude  larger  than  the  pure  delay. 

The  number  of  nodes  needed  for  a  good  approximation  appears  to  be  very  small  [55].  In 
particular,  models  with  two  nodes  might  be  adequate,  although  inputs  with  high  bandwidth 
clearly  necessitate  a  few  more. 
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2.8  Summary 


This  chapter  began  with  the  complete  equations  of  motion  for  a  cable  in  water,  and  detailed  a 
number  of  valid  simplifications.  These  pertain  to  the  neglect  of  torsional,  rota*.  ^  >nal  inertia, 
and  elongation  effects,  and  to  the  specific  forms  of  the  added  mass  and  quadratic  drag 
terms.  In  addition,  geometric  coupling  was  largely  eliminated.  Finally,  since  a  very  low 
order  finite-difference  approximation  is  advocated,  the  continuous  nature  of  the  cable  is 
lost  as  well.  The  result  is  a  simple  model  structure  suitable  for  control  system  design  and 
analysis. 

The  lateral  frequency  response  of  a  deeply-towed  cable  system  is  strongly  dependent 
on  the  amplitude  of  excitation,  as  shown  by  both  the  harmonic  balance  and  perturbation 
studies  of  this  chapter.  The  general  rule  of  thumb  coming  out  of  the  analysis  is  that  the  gain 
of  the  plant  is  inversely  related  to  the  square  root  of  the  excitation  frequency.  In  addition, 
the  analyses  reveal  a  deterioration  of  resonant  peaks  with  growing  excitation  amplitude, 
and  an  arbitrarily  large  loss  of  phase,  even  well  below  the  first  natural  mode.  These  results 
are  suggestive  of  the  need  for  nonlinear  control  strategies. 

It  was  then  shovm,  however,  that  control  spillover  effects  are  likely  to  be  damped  out 
very  quickly  in  practice,  due  to  the  large  drag  terms.  Also,  closed-loop  limit-cycling  brought 
about  by  application  of  a  linear  feedback  law  was  predicted  to  be  of  quite  small  sunplitude, 
even  for  a  poorly-designed  compensator.  These  results  bode  well  for  linear  control,  in  a 
very  limited  sense. 
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Chapter  3 


Open-Loop  Maneuvers 


3.1  Introduction 

Studying  maneuvers  in  which  the  vessel  motions  are  “preshaped”  is  important  as  a  first  step 
in  positioning  the  towed  imderwater  system  more  accurately.  Preshaping  is  an  open-loop 
idea  and  quite  distinct  from  feedback  control,  which  is  discussed  in  Chapter  4.  Preshaping 
accomplishes  our  first  objective,  which  is  to  speed  up  the  dynamic  response  of  the  towed 
system,  by  anticipating  in  a  constructive  way  how  the  plant  will  respond  to  excitation. 
Because  it  is  an  open-loop  procedure,  it  h2is  no  instabilities,  and  is  therefore  “safe”  from 
the  point  of  view  of  the  helm. 

We  discuss  two  aspects  of  preshaping  vessel  trajectories  in  this  chapter  The  first  involves 
the  so-called  point-to-point  maneuver:  the  object  is  to  move  the  lower  end  of  the  cable  or 
riser,  at  rest,  to  a  new  location  and  achieve  the  rest  configriration  again  in  miTiiTmiTo  time. 
The  second  approach  involves  specified  trajectories  for  the  towfish,  and  more  mathematics  to 
invert  the  cable  d3rnamics.  In  this  chapter,  we  outline  the  techniques,  and  present  verifying 
experimental  data  from  a  full-scale  test. 

3.2  Point-to-Point  Moves 

The  point-to-point  move  is  a  common  one  for  many  operations,  and  has  been  considered 
numerically  by  Paul  and  Soler  [80].  Neglecting  the  presence  of  currents  for  the  moment,  we 
assume  that  the  system  begins  at  rest;  that  is,  the  tether  is  hanging  vertic2dly  due  to  the 
weight  of  the  system.  The  goal  is  to  move  the  towfish  horizontally  to  a  new  location  and  to 
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Fig\ire  3-1:  Slow  towfish  response  to  a  point-to-point  vessel  maneuver. 


bring  it  to  rest,  by  moving  the  surface  vessel  appropriately.  As  was  seen  in  Chapter  2,  the 
cable  usually  acts  like  a  low-pass  nonlinear  filter  in  transmitting  ship  motions  to  towfish 
motions.  This  property  is  illustrated  in  another  way  in  Figme  3-1,  which  shows  a  simulated 
response  for  the  case  where  the  vessel  is  moved  directly  to  a  new  location.  The  towfish 
settling  time  to  10  meters  for  this  example  (2000  meters  of  6.4-millimeter  steel  wire,  740 
Newtons  static  tension  at  the  bottom)  is  approximately  40  minutes. 

It  is  intmtively  clear  that  to  reduce  the  settling  time  of  the  fish,  we  should  alter  the 
vessel  trajectory  to  incorporate  the  cable  response  more  vigorously.  There  are  many  ways 
of  posing  this  problem  mathematically;  we  consider  the  teisk  of  minimizing  the  time  integral 
of  towfish  position  error.  More  formally,  let  the  function  g  map  the  vessel  trajectory  Zv(t) 
to  the  towfish  trajectory  Then  we  want  to  satisfy 

fif 

*„(<)  =  arg  min  /  |  I  (3-1) 

Jto 
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subject  to  various  constraints  on  the  vessel  motions.  Note  that  Zt,ile««rc<l  is  ^  step  function 
for  the  point-to-point  maneuver. 

3.2.1  Approximate  Solution 

Given  that  g  in  Equation  3.1  is  a  nonlinear  operator  with  memory,  there  appears  to  be 
uo  general  closed-form  solution  to  the  optimization  problem  stated  above.  However,  it 
can  be  argued  that  the  best  vessel  trajectory  moves  the  whole  cable  first,  and  the  towfish 
indirectly.  Thus,  we  consider  an  in'tial  vessel  overshoot  to  “stretch  out”  the  cable  over  the 
desired  location,  and  then  a  return  to  that  location.  During  this  return,  the  vehicle  falls  into 
place  over  the  target.  This  is  the  basic  approach  taken  by  Paul  and  Soler,  who  minimized 
the  settling  time  to  a  given  radius.  Since  the  overshoots  of  this  thesis  and  of  their  paper 
have  the  same  shape,  it  is  clear  that  the  two  minimization  problems  are  very  similar;  the 
settling  time  minimization  of  of  Paul  and  Soler’s  paper  requires  an  extra  parameter,  namely 
the  settling  radius. 

Our  general  overshoot  maneuver  may  be  parameterized  by  four  variables,  shown  in 
Figure  3-2: 

1-  Omax'-  the  absolute  value  of  the  maximum  acceleration  of  the  vessel 

2.  V/:  the  absolute  value  of  the  maximum  forward  speed  of  the  vessel 

3.  Vr‘  the  absolute  value  of  the  maximum  reverse  speed  of  the  vessel 

4.  treverje^  the  time  elapsed  at  which  the  vessel  changes  from  forward  to  reverse  move¬ 
ment 

The  acceleration  and  speeds  are  generally  governed  by  operational  considerations;  e.g.,  how 
hard  to  work  the  DP  system,  or  what  cable  angles  are  acceptable.  For  many  operations,  the 
maximum  speeds  are  only  several  knots.  The  distinction  between  maximum  forward  and 
reverse  speeds  is  made  because  the  cable  may  be  drawn  close  to  the  vessel  hull  when  backing 
down.  Thus,  we  have  effectively  reduced  the  number  of  unknowns  in  the  trajectory  x„{t)  to 
one:  the  reversal  time  irevtrtt  *  The  problem  is  reduced  to  a  single  parameter  optimization 
of  the  cost  integral  of  Equation  3.1. 

Paul  and  Soler  showed  that  the  settling  time  for  the  towfish  is  a  convex  function  of 
the  overshoot  distance;  this  implies  that  lengthy  search  algorithms  are  not  necessary  to 
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®max 


Figure  3-2:  Assumed  optimal  vessel  trajectory. 


find  trtversf  Instead,  we  may  simply  start  from  the  zero  overshoot  trajectory,  and  increase 
Utverie  Until  the  computed  cost  integral  has  reached  its  miniTniiTn.  The  compactness  of 
our  finite-difference  model  plays  a  crucial  role  in  malting  this  procedure  computationally 
feasible. 

3.2.2  Full-Scale  Experimental  Setup 

The  full-scale  data  to  be  presented  were  obtained  during  a  cruise  at  the  Juan  de  Fuca  Ridge 
in  July  1991,  aboard  R/V  Laney  Cbouest.  This  ridge  is  in  the  Pacific  Ocean,  approximately 
300  kilometers  west  of  Washington  state,  and  has  been  the  site  of  numerous  recent  geological 
studies.  The  deployment  was  designed  as  a  fairly  close  approximation  of  the  ARGO/JASON 
system,  and  in  fact  was  set  up  just  prior  to  ARGO/JASON  operations;  it  is  schematized 
in  Figure  3-3  and  Table  3.1.  The  vessel  was  navigated  with  a  dynamic  positioning  system, 
using  fixes  from  the  Global  Positioning  Satellite  system  and  a  long-baseline  acoustic  net. 
The  acoustic  frequencies  were  between  five  and  fifteen  kilohertz,  and  the  net  had  an  average 
side  length  of  approximately  five  kilometers.  An  acoustic  transceiver  was  placed  at  the  end 
of  2000  meters  of  hydrographic  wire,  along  with  a  steel  clump  weight,  and  the  motions  of 
this  clump  weight  were  tracked  with  the  same  long-baseline  net.  The  vessel  was  headed  into 
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Figure  3-3:  Schematic  view  of  full-scale  deployment. 

the  sea  during  each  of  the  runs,  and  the  weight  was  deployed  over  the  side  amidships.  There 
were  no  significant  currents  in  any  of  the  experiments  conducted,  based  on  the  observed 
offsets. 

The  required  vessel  trajectories  were  executed  with  the  existing  dynamic  positioning 
system.  Since  this  system  did  not  have  the  capability  to  follow  complex  trajectories,  we 
operated  in  regulation  mode  and  had  the  setpoint  foUow  the  trajectories  specified  by  our 
algorithms.  Our  model  of  the  experimental  plant  was  initially  estimated  from  its  physical 
parameters  (mass,  weight,  cable  diameter,  etc.;  see  Section  2.7),  before  the  crtiise.  This 
model  was  then  updated  at  sea  using  preliminary  test  data,  and  the  experiments  which 
follow  were  based  on  the  new  model,  which  had  four  nodes.  Finally,  the  data  which  follow 


70 


Actuator 

R/V  Laney  Chouest  (2400  tons) 

Maximum  towing  speed 

1  m/s 

Vessel  tracking 

GPS,  acoustic  (<r  ^  5  m) 

Vessel  DP  bandwidth 

.02  rad/s 

Cable  material 

6.4  mm  hydrographic  wire 

Cable  effective  mass 

.181  kg/m 

Cable  weight /length  in  water 

1.14  N/m 

Cable  Cd 

1.4 

Depth 

2000  m 

Towfish 

steel  plates  (740  N) 

Towfish  effective  mass 

100  kg 

Towfish  frontal  area 

340  cm* 

Towfish  Cd 

2 

Towfish  tracking 

acoustic  (<7  c::  5  m) 

Pure  delay 

:::  20  s 

Table  3.1:  Physical  parameters  for  full-scale  experiments. 


have  not  been  post-processed  in  any  way,  except  for  the  manual  removal  of  several  “fliers.” 

3.2.3  Overshoot:  Full*Scale  Experimental  Results 

As  described  above,  the  goal  of  the  experiment  was  to  show  that  point-to-point  relocation 
of  the  towfish  can  be  made  more  efficient  by  using  the  calculated  overshoot  idea.  Our 
simxilations,  and  those  of  Paul  and  Soler  [80],  show  that  reductions  in  the  towhsh  settling 
time  (to  a  specified  radius)  can  approach  twenty-five  percent  or  more,  depending  on  the 
physical  system  and  the  parameters  of  the  move.  We  give  two  experimental  illustrations  of 
this  reduction  in  settling  time. 

Figure  3-4  shows  an  overlay  of  results  from  two  tests;  one  test  consisted  of  a  flat  vessel 
move  at  0.25  meters /second,  and  the  other  had  an  overshoot  of  approximately  150  meters, 
at  0.25  meters/second.  The  relocation  distance  was  300  meters  for  both  runs.  The  settling 
time  to  10  meters  for  the  flat  case  is  approximately  45  minutes,  while  for  the  overshoot, 
it  is  approximately  40  minutes  (this  takes  into  account  the  fact  that  our  vessel  trajectory 
settled  at  310  meters). 

The  second  case  involved  a  much  more  vigorous  overshoot  for  the  vessel;  see  Figxire  3-5. 
The  relocation  here  was  100  meters,  with  maximum  velocities  of  0.5  meters/second.  For  this 


71 


Figure  3-4:  Comparison  between  overshoot  and  flat  responses  for  flrst  full-scale  move. 

run,  we  let  the  vessel  stay  relatively  stationary  for  part  of  the  “stretch”  period,  instead  of 
making  a  single  sharp  reversal  at  trcver«e-  The  flat  run  has  a  10-meter  settling  time  of  about 
22  minutes,  although  there  are  apparent  transients  at  the  start  and  at  the  end  of  the  data 
set.  The  overshoot  leads  to  a  different  response,  characterized  by  faster  towflsh  motion  as 
before,  and  a  much  sharper  stop.  Note  that  through  modeling  and  vessel  trajectory  errors, 
the  towflsh  itself  has  an  overshoot  of  approximately  10  meters.  The  10-meter  settling  time 
for  the  overshoot  is  approximately  17  minutes;  the  difference  between  the  flat  and  overshoot 
responses  is  more  significant  for  this  nm. 

Another  way  to  see  the  effects  of  the  overshoot  is  to  look  at  a  side  view  of  the  towflsh 
motion;  Figme  3-6  gives  such  a  view  for  the  nms  of  Figure  3-5).  We  see  that  the  towflsh 
trajectories  for  the  flat  and  overshoot  cases  are  nearly  identical  early  in  the  move,  with  high 
frequency  vertical  oscillations  caused  by  heaving  of  the  surface  vessel.  Later,  while  the  flat 
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Figure  3-5:  Comparison  between  overshoot  and  flat  responses  for  second  full-scale  move. 
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Figure  3-6:  Side  views  of  overshoot  and  fiat  towfish  trajectory. 

run  requires  some  time  for  the  towfish  to  drift  slowly  down  into  the  new  rest  position,  the 
overshoot  leads  to  a  towfish  trajectory  which  is  better  described  as  dropping.  This  behavior 
is  consistent  with  the  notion  that  the  overshoot  move  is  trying  to  stretch  the  cable  out  over 
the  desired  rest  position,  and  then  let  it  drop  into  place. 

As  a  more  practical  use  for  the  overshoot  concept,  we  now  consider  the  construction  of 
patterns  (in  the  horizontal  plane)  using  cascaded  overshoot  maneuvers  in  difierent  direc¬ 
tions.  The  primary  motivation  for  this  extension  is  grid  searches  in  which  the  towfish,  often 
gathering  acoustic  data,  must  move  in  straight  lines  in  a  grid.  The  traditional  procedure 
involves  a  lengthy  turnaround  maneuver  at  the  end  of  each  tracldine,  which  aUows  the  ca¬ 
ble  to  become  aligned  in  one  vertical  plane  before  starting  the  next  trackline.  With  our 
overshoot  procedure,  sharp  comers  in  the  towfish  trajectory  can  be  achieved,  eliminating 
the  necessity  of  the  standard  turnaround  maneuvers.  See  Figure  3-7. 
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Figure  3-7:  Cascaded  overshoots  for  precision  grid-following 
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Figure  3-8:  Bird’s-eye  view  of  sharp  comers  executed  by  the  towhsh. 

Three  overshoot  maneuvers  were  pieced  together  in  immediate  succession  to  obtain  the 
full-scale  results  of  Figure  3-8,  which  gives  a  bird’s-eye  view  of  the  vessel  and  resulting 
towfish  trajectories.  The  first  leg  moved  the  towfish  550  meters  (at  0.5  meters/second)  to 
the  southeast,  the  second  leg  moved  it  100  meters  (at  0.25  meters/second)  to  the  southwest, 
and  the  third  leg  moved  it  200  meters  (at  0.5  meters/second)  to  the  northwest.  These 
projections  are  given  in  Figure  3-9.  All  three  of  the  overshoots  were  effective,  especially  the 
two  northwest-southeast  moves. 

3.2.4  Discussion 

In  the  first  two  in-plane  runs,  the  approximate  reductions  in  settling  time  to  ten  meters  are 
approximately  11%  and  23%  of  the  flat  maneuver  settling  time;  these  numbers  are  consistent 
with  those  predicted  by  simulations.  The  overshoot  maneuver  appears  to  be  most  useful  in 
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meters 


Projection  on  300Klegree  azimuth 


Figure  3-9:  Projections  along  perpendicular  azimuths  for  three  overshoot  moves. 


77 


moves  of  short  net  distance  and  high  vessel  speed,  also  consistent  with  previous  results. 

The  most  dramatic  effects  of  the  overshoot  maneuver  are  found  in  the  cascaded  nm  of 
Figure  3-8,  however.  Since  the  emphasis  is  on  turning  sharp  comers  with  the  towfish,  it 
is  interesting  to  compare  the  experimental  results  to  the  response  of  a  model  without  the 
overshoots  (in  lieu  of  experimental  results  with  flat  vessel  trajectories).  To  this  end,  we 
used  a  model  with  coupled  drag  terms,  and  drove  it  with  the  recorded  vessel  trajectory  for 
this  run,  but  with  the  overshoots  truncated.  A  zoom  view  of  the  experimental  comers  and 
the  sim\ilation  comers  is  given  in  Figure  3-10.  We  see  that  the  simulation  with  flat  vessel 
motion  is  missing  the  sharp  comers  found  in  the  experimental  data,  to  the  order  of  fifty 
meters  error,  and  yet  the  time  to  execute  the  total  maneuver  is  identical  in  both  cases. 
The  implication  is  that  following  a  trackline  with  corners  cannot  be  achieved  as  precisely 
without  the  overshoots,  for  the  same  time  period  and  the  same  vessel  speeds.  Indeed,  the 
deterioration  of  this  northwest-southeast  trackline  is  a  good  example  of  why  the  usual  long 
turnaround  has  traditionally  been  used. 

An  important  property  of  the  overshoot  response  is  that  the  towfish  goes  through  sub¬ 
stantial  altitude  changes;  during  the  vessel  “return”  period,  the  rate  of  altitude  loss  may 
be  high.  For  the  experiments  reported  here,  the  maximum  vertical  velocity  of  the  towfish 
was  under  0.2  meters/second.  Many  winches  can  keep  up  with  velocities  in  this  range,  and 
it  is  expected  that  constant  altitude  during  these  overshoot  moves  can  be  maintmned,  if 
necessary,  by  actively  using  the  winch.  One  must  be  aware,  however,  that  the  vessel  tra¬ 
jectories  of  this  section  do  not  involve  horizontal  velocities  greater  than  0.5  meters /second; 
the  altitude  problem  will  surely  be  exacerbated  in  operations  at  higher  speeds. 


3.3  Open-Loop  Trajectory  Following 

We  consider  now  the  task  of  input  preshaping  for  more  general  trajectory  following.  As¬ 
suming  that  the  operator  has  specified  the  trajectory  of  the  towed  body  as  a  time  series,  the 
problem  is  to  specify  the  required  ship  motions  to  bring  about  this  response.  The  special 
case  of  a  time-optimal  step  change  in  position  of  the  lower  endpoint  (as  above)  is  a  subclass 
of  the  trajectory-following  problem. 
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Figure  3-10:  Comparison  of  corners  achievable  with  (top)  and  without  (bottom)  overshoots. 
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3.3.1  Background 

The  concept  of  input  preshaping  is  an  intuitive  one;  we  seek  to  “undo”  the  dynamics  of 
the  physical  system,  calling  what  is  an  output  from  the  real  system  (towfish  motion)  an 
input  to  the  algorithm,  and  vice  versa.  The  alternative  to  preshaping  is  feedback  applied 
to  a  moving  reference  point.  This  technique  will  be  considered  experimentally  in  Chapter 
5,  but  it  is  clear  that  robustness  considerations  are  likely  to  cause  sluggish  performance  in 
comparison  to  a  full  inversion  procedure. 

A  novel  approach  has  recently  been  applied  to  the  feedforward  component  of  tracking 
control  for  flexible  robots  by  Bayo  et  al.  [9],  [11],  [12],  [10].  In  short,  they  use  a  finite- 
element  discretization  of  the  link,  and  an  inversion  in  the  frequency  domain.  Given  the 
desired  acceleration  profile  of  the  endpoint,  and  if  the  arm  follows  linear  dynamics,  a  single 
computation  gives  the  necessary  hub  torque  profile.  One  notable  feature  of  this  approach  is 
that,  because  the  inversion  occurs  in  the  frequency  domain,  the  hub  torque  command  can 
precede  the  desired  endpoint  motions  so  as  to  accoimt  for  the  flexure  in  the  arm,  and  the 
piire  time  delay.  In  the  case  of  nonlinear  terms  in  multilink  arms  (Coriolis  cmd  centrifugal 
torques),  these  researchers  have  implemented  a  simple  iterative  scheme  which  apparently 
converges  to  the  correct  results. 

We  begin  with  this  formulation  and  modify  it  to  fit  otir  system.  The  towing  problem 
differs  in  that  pointwise  nonlineaxities  must  be  allowed,  since  the  cable  follows  a  quadratic 
drag  law.  In  addition,  a  flnite-difference  scheme  is  used  to  model  the  distributed  system, 
instead  of  the  finite-element  method;  these  models  are  described  in  Section  2.7.  The  only 
other  change  is  that  we  initially  neglect  bending  stiffness  for  the  simple  reason  that  it  is 
small  for  our  applications. 

In  the  two  sections  that  follow,  we  outline  the  iterative  method  of  inversion  for  the 
in-plane  problem,  and  then  give  a  proof  for  the  convergence  of  the  nonlinear  scheme. 

3.3.2  The  Iterative  Method 

In  the  following,  the  qualifier  “(t)”  indicates  a  time-domain  qu£uitity,  and  “(u;)”  indicates 
the  frequency  domain  equivalent.  Lk  is  a  column  of  zeros  with  a  one  in  the  l;-th  row,  and 
I-k  is  the  identity  matrix  with  a  zero  replacing  the  one  in  the  k-th  row.  The  sizes  of  tj,  and 
J-k  will  be  clear  from  context. 
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The  generalized  n-node  system  is  modeled  as  in  Section  2.7,  but  we  cast  it  in  a  different 
form: 


+  £*'(<)  +  Kx(t)  =  Gu(t)  -  Nit,  *(t),  x\t)),  (3.2) 

where  M,  B,  and  K  are  the  linear  mass,  damping,  and  stiffness  matrices,  respectively. 
M  and  B  are  diagonal,  and  K  is  tridiagonal  for  the  cable  system  without  bending  stiff¬ 
ness.  The  vector  z(t)  is  the  length  n  -f-  1  position  vector  including  the  vessel;  x(t)  := 
[z„(t)  qn-i  •  *  ■  9i]^.  This  notation  makes  qi  the  position  of  the  towfish,  consistent 
with  the  discretization  scheme  of  Chapter  2.  u(t)  is  the  input  (in  units  of  force)  at  the  ves¬ 
sel,  so  that  G  =  L\.  JV(t, x(t),z'(t))  is  the  nonlinear  force  vector  to  be  specified.  Without 
loss  of  generality,  we  can  write  Equation  3.2  in  the  frequency  domain  to  give 

+  iwB  -I-  A’]z(u>)  =  Guiu))  -  Af{u).  (3.3) 

Here,  z(a>)  is  the  Fourier  trauisform  of  the  position  vector  z(t).  A/'(a»)  can  be  thought  of  as 
the  Fourier  transform  of  the  time-domain  nonlinear  force  Nit).  Renaming  the  bracketed 
expression  H{u))  and  setting  Fiu)  :=  J7~^(a>),  one  finds  that  the  linear  solution  comes  from 


z(u>)  =  f  (u;)^^^),  or 

(3-4) 

When  O’  =  0,  a  diagonal  matrix  with  small  elements  can  be  added  to  K  to  make  ■H(O) 
invertible.  In  addition,  Herbert  and  Jones  [54]  noticed  that  when  Hiu)  is  tridiagonal, 
Fn+i,ii^)  *  particularly  simple  cofactor  form,  which  reduces  the  inversion  of  Hiu>)  to 
essentially  one  determinant  calculation.  This  fact  would  be  useful  for  a  very  high-order 
system. 

Thus,  the  procedure  in  the  linear  case  is  to  first  specify  91  (t)  as  a  finite  time  series 
of  length  p,  and  convert  it  to  the  frequency  domain.  The  above  equation  then  has  to  be 
solved  for  each  distinct  frequency  in  the  domain  of  interest  for  the  physical  system.  Once 
the  transformed  input  has  been  solved  for  the  relevant  frequency  range,  one  inverse  Fourier 
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Transform  gives  the  required  input  for  the  physical  system.^ 

In  the  lumped-parameter  framework,  we  propose  that  the  nonlinear  term  may  be  itera¬ 
tively  solved  as  follows: 

1.  Derive  an  expression  for  u((i;}  involving  the  linear  part  found  above,  and  the  nonlinear 
part  (initially  zero).  There  are  at  least  several  ways  to  do  this — the  form  we 

use  is  as  follows: 


«i+i(w)  =  (3.5) 

where  the  k  subscript  denotes  the  iteration  index.  This  formula  is  the  simplest  to 
derive,  and  has  a  nonlinear  part  which  is  additive  with  the  original  linear  component. 

2.  Evaluate  the  pointwise  position  transforms  based  on  this  new  u(a>).  Note  that  z„.^i(a>) 
does  not  need  to  be  reevaluated  since  it  is  just  qi .  We  have 


**+1(4^)  =  ‘n+l9l(‘^)  +  •^-(n+l)'f’(‘*')(^«fc+l(w)  -  Mt(u>))-  (3.6) 

3.  Calculate  the  new  value  of  the  nonlinearity  based  on  **^.1  (u>);  all  time  derive  '; '  vt  s  ad 
integrals  are  related  by  multiples  of  the  factor  iu>.  We  bundle  these  operations,  in 
addition  to  the  fundamental  nonlinearity  term,  into  the  operator 


A/J,+i(a;)  =  ^(u»,*fc+i(fa>i),Zfc+i(u>2),...,z*+i(u;p)).  (3.7) 


As  written,  the  nonlinearity  may  involve  more  than  one  of  the  p  frequency  components. 
For  example,  a  square  nonlinearity  makes  <f>  a  convolution  operator  in  the  frequency 
domain. 


similar  relation,  which  avoids  the  finite-difference  approximation,  can  be  derived  from  the  linear 
wave  equation,  with  appropriate  boundary  conditions.  It  takes  the  form  of  an  irrational  transfer  function 
from  the  input  u  to  the  output  91;  to  solve  the  inverse  dynamics,  one  merely  has  to  invert  the  transfer 
function.  However,  this  approach  is  somewhat  cumbersome  with  depth- varying  tension  and  is  impractical 
when  nonlinearities  are  introduced. 
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Figtire  3*11:  Iterative  loop  as  a  discrete- time  system. 


3.3.3  Stability  Proof  for  the  Iterative  Scheme 

The  procedure  of  the  previous  section  can  be  seen  as  a  discrete  feedback  loop  involving  a 
new  and  very  large  (n  +  l)p  x  1  state  vector  X  :=  [a:(u>i)^  This 

loop  is  drawn  in  Figure  3-11,  and  is  described  by 

=  Fxn+i  +fii{Xk),  (3,8) 

where  x«  :=  [®»(wi)  a;i(w2)  •••  F  is  an  (n-|- l)p x p  constant  matrix,  and /3  is  an 

(n  +  l)p  X  (n  -f-  l)p  constant  matrix.  We  have  also  set 

^(wi,*(a»i a:(a;p)) 

i{X):=  : 

^(wp,x(a;i),--.,a:(u;p)) 

The  nonlinear  operator  i  carries  into  itself.  It  is  noteworthy  that  we  have  totally 

eliminated  u  from  the  loop;  it  is  only  ii  :=  *„  that  is  of  interest. 

In  order  to  show  loop  stability,  we  consider  the  transition  matrix  /3,  and  then  its  interaction 
with  the  nonlinearity  i.  To  begin,  ^  is  the  block-diagonal  of  p  frequency-specific  parts 

0 

For  this  system,  the  matrix  is  completely  defined  by  and  some  trivial  vectors 
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Figure  3-12:  /3  as  a  matrix  of  two-determinants  of  F. 

as  follows: 

/3'(a,,)  =  j  ,  (3.9) 

and  on  carrying  through  the  m\iltiplications,  we  find  that  has  a  first  column  and  last 

row  of  zeros.  Dropping  the  frequency  argument  for  clarity,  the  other  terms  are  given  by: 

-  FiJ  for  i  =  1, n  and  i  =  2, n  4- 1.  (3.10) 

•^n+1,1 

We  are  particularly  interested  in  the  elements  of  )  in  the  lower  triangle;  these  terms 
are  multiples  of  specific  two-determinants  of  F(u>i)  which  are  cornered  in  the  lower  left  of 
F(u>;).  For  this  reason,  these  elements  of  will  be  called  “corner-2-determinants”. 

Their  graphical  interpretation  is  given  in  Figure  3-12. 

The  following  proposition  applies  to  a  single  H(ui),  and  hence  the  argument  Ui  is  omitted 
again  for  clarity: 

Proposition:  Let  if  be  an  invertible  tridiagonal  matrix  of  size  (n-l- 1)  x  (n-f- 1).  The  comer- 
two-determinants  described  above,  which  involve  the  lower  left  comer  element  of  H~'^  and 
any  other  element  which  lies  in  the  lower  triangular  portion  of  are  identicaUy  zero. 
More  formally, 

Fi,\Fn+\,j  —  =  0  for  t  =  2, •  •  -  jH  and  7  =  2,  •  •  •,  t.  (3-11) 
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Proof:  The  algebraic  proof  proceeds  in  three  parts:  we  show  successively  that  the  comer- 
2-deterniinants  of  the  1)  first  and  second  columns,  2)  of  the  first  and  third  columns,  and  3) 
inductively  of  the  first  and  j-th  columns,  are  zero. 


1.  The  definition  of  the  inverse  matrix  is  FH  =  I.  This  gives 

+  -ft.r-ffj.i  =  0,  t  =  2,  •  •  • ,  n  +  1.  (3*12) 

One  obtains  the  desired  result  immediately  for  the  second  column  of  /3',  by  considering 
the  i-th  and  the  n  +  1-th  rows; 

=  0  for  1  =  2,  •  •  • ,  Tl.  (3.13) 

2.  The  identity  gives,  for  the  remaining  lower  triangle  elements  of  F: 

4-  Fi,jHjj.i  =  0  for  >  =  3,  •  •  • ,  n -j- 1,  t  =  ;,  •  •  • , n -f  1. 

(3.14) 

Consider  the  case  j  —  3;  we  have 

+ -fi.s-ffs.r  =  0  for  t  =  3,  •  •  • ,  n,  (3.15) 

■fn+l,l-ffl,2  +  •fn+l,2-ff2,2  +  Fn+l,3-ff3,2  =  0.  (3.16) 

We  substitute  the  value  of  Fi,2  from  Equation  3.13  into  Equation  3.15,  multiply  Equa- 

—  IT- 

tion  3.16  by  1  and  add  the  results.  This  gives 

-f’i,l-fn+l,3  —  -fi.S-fn+l.l  =  0  for  t  =  3,  •  •  • ,  Tl.  (3.17) 

3.  Now  we  generalize  to  the  case  for  arbitrary  j,  but  we  must  start  with  j  =  4,  and 

incrementally  move  to  j  =  n  -f-  1  in  Equation  3.14.  As  in  Part  2),  we  can  write  a  pair 

of  equations,  one  for  the  t-th  row,  and  one  for  the  n  +  1-th  row: 
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0,  i  =  j, 
0. 


,n,(3.18) 

(3.19) 


In  general,  the  computations  for  the  previous  two  columns  provide  the  ‘‘bridge”  back 

to  the  first  column.  For  example,  in  the  case  j  =  4,  one  substitutes  the  value  of 

Fi,3  from  Equation  3.17  and  Fi,2  from  Equation  3.13  into  Equation  3.18,  multiplies 
“F' 

Equation  3.19  by  as  before,  and  adds  the  results.  The  conditions  on  the  previous 
two  columns  are 


■fi.ifn+ij-i  -  =  0  for  t  =  j  -  1,  •  •  • ,  n,  (3.20) 

Bi.iFn+ij-i  -  =  0  for  i  =  J  -  2,  •  •  • , Ti.  (3.21) 

These  operations  are  repeated  until  y  =  n  is  reached;  in  this  way,  the  proposition  is 
proved.  □  □  □ 

Each  P'{u!i),  and  thus  /3,  has  eigenvalues  which  are  zero,  since  the  lower  triangular 
elements  are  zero.  As  will  be  shown,  however,  /3  has  a  particularly  interesting  structure 
which  protects  the  discrete-time  system  d}rnamics  from  potentially  destabilizing  nonlinear 
effects  and  guarantees  convergence  in  n  iterations.  We  need  the  following  assumption:  the 
nonlinearity  #  may  span  across  different  frequency  components,  but  does  not  span  &om 
node  to  node  in  the  physical  system.  Mathematically,  this  means 

dii 

=  0  unless  t  =  y  (n  -|-  l)k,  k  some  integer.  (3.22) 

0^  j 

This  is  a  sufficient  condition  for  the  final  thrust  of  the  section:  it  specifies  that  the  nonlinear 
effects  depend  on  physical  quantities  which  are  local.  The  fluid  drag  nonlinearity  fits  this 
assumption  neatly,  but  a  Coriolis  torque  on  a  flexible  manipulator,  for  example,  does  not. 

With  the  above  assumption  and  the  known  structure  of  /3,  it  is  apparent  that  the  discrete 
system  of  Figure  3-11  will  converge  in  n  iterations.  To  show  this,  we  make  the  following 
arguments.  First,  yd  represents  an  order-n  deadbeat  system  matrix,  since  its  only  nonzero 
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elements  are  in  the  diagonal  blocks  /3'(a;,-),  which  are  themselves  order-n  deadbeat  system 
matrices  (order-n  because  ^  does  not  affect  Xn+i)*  Thus,  if  i  were  set  to  the  identity,  an 
additional  p  states  in  X  would  become  invariant  on  each  iteration,  giving  total  convergence 
of  ^  in  n  iterations. 

In  the  presence  of  a  local  nonlinearity  $,  the  nature  of  this  convergence  is  tmchanged. 
The  reason  is  that  each  element  of  #(X)  depends  on  only  one  Xit  consistent  with  Equa¬ 
tion  3.22.  This  dependence  is  so  arranged  that  the  order-n  deadbeat  nature  of  is  fully 
exploited. 

Example:  Suppose  n  -|- 1  =  3  and  p  =  2;  then 


0  fii,3  0  0  0 

0  0  02,3  0  0  0 

0  0  0  0  0  0 

0  0  0  0  04^s  04,6 

0  0  0  0  0  /Ss.e 

0  0  0  0  0  0 


The  lower  endpoint  motions  are  fixed,  so  also  fixed  are  X3>  *3  and  *6-  We  have  Xi  =  rx3, 
and  let  x^  be  the  original  value  of  Xi  and  Zi  be  the  origined  value  of  z^.  With  some  abuse 
of  notation,  the  evolution  of  the  state  X  is  as  follows: 


^2 


0l,2^2iX2)  +  /3i,3^3(X3)  +  *1 
02,3^3(X3)  +  *2 
*3 

^4,5^s(X2)  +  134,6$6(X3)  +  *4 
^5.6^6(X3)  +  *5 
*6 
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^1,J^2(/32,3^3(X3)  +  *2)^6,«^6(X3)  +  *$)  +  ^1,3^3(X3)  +  *1 
/32,3*3(X3)  +  *2 
£3 

X3=  ,  and 

^4,5^5(^2,3^3(X3)  +  *2»^6.6^6(Xs)  +  *b)  +  04,6^e(X3)  +  *4 

fis.e^eiXs)  +  *5 

*6 

01, 2^2(02, 3^3(X3)  +  *2,^5.6^«(X3)  +  *5)  +  /3i.3#3(X3)  +  *1 
;32,3*3(X3)  +  *2 
*3 

0A,3^s{02,3^3{X3)  +  i2,/36.6$6(X3)  +  *5)  +  04,6^3{X3)  +  ®4 
^S.6*6(X3)  +  *5 
*6 

One  sees  that  X3  is  indeed  unaffected  by  the  loop.  After  the  first  iteration,  it  is  clear  that 
X2  is  fixed,  and  after  the  second,  that  xi  is  fixed. 

The  iteration  scheme  generates  the  complete  motion  of  successive  points  along  the  ca¬ 
ble,  one  at  a  time.  This  is  precisely  how  one  would  solve  this  problem  with  backward  time 
integration,  one  node  at  a  time.  However,  it  has  been  shown  by  Serna  and  Bayo  [88]  that 
straightforward  backwards  integration  cannot  work  in  general,  due  to  numerical  instabili¬ 
ties.  A  time-domain  computation  of  the  impulse  response  integral  has  been  suggested  as 
a  panacea,  although  such  an  approach  is  valid  only  for  linear  systems.  We  do  not  pmsue 
time-domain  computations  further  here. 

3.3.4  Applications 

Many  physicaJ  systems  fit  into  the  tridiagonal  system  matrix  form,  and  lateral  spring-mass- 
damper  approximations  of  the  wave  equation  are  clearly  suitable.  The  assumption  3.22 
allows  a  very  large  range  of  nonlinearities  to  act  upon  such  a  system,  as  long  as  they  act 
locally,  and  the  an2ilysis  indicates  that  the  rate  of  convergence  is  independent  of  the  nature 
of  this  nonlinearity.  The  ciirrent  section  presents  four  examples  (simulations  only)  of  the 
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usefulness  of  this  procedure. 

An  important  note  is  that,  for  ocean  applications,  we  require  the  Fourier  series  to  be 
truncated.  This  is  in  contrast  to  Bayo*s  results,  in  which  torque  commands  with  very 
high  bandwidth  can  be  applied.  An  immediate  reason  to  truncate  is  that  a  reduction  in 
computation  time  by  several  orders  of  magnitude  can  be  gained,  since  there  are  several 
inversions  of  for  each  frequency  retained.  Another  consideration  is  that  when  the  full 
series  is  kept,  the  algorithm  presented  here  may  give  discontinuities  or  very  high  velocities 
in  Zv,  even  if  qi  is  smoothed.  At  present,  we  have  no  mechanism  in  place  that  prefers 
smooth  Zv!  obviously,  one  could  filter  the  resultant  x„  after  the  algorithm  is  completed. 

Extensive  experience  with  the  algorithm  indicates  that  stability  of  the  algorithm  is  an 
issue  only  when  the  number  of  terms  in  the  Fourier  series  is  inconsistent  with  the  order  of 
the  plant.  For  example,  it  is  physically  impossible  to  obtain  fast  transients  at  the  towfish 
if  the  cable  model  is  second-order.  The  higher  the  order  of  the  plant,  the  sharper  are 
the  transients  we  are  able  to  achieve,  using  more  frequency  components.  In  addition,  the 
algorithm  appears  to  fail  if  any  part  of  the  cable  takes  on  this  imbalance.  The  typical  failure 
scenario  is  as  follows*,  a  smooth  towfish  trajectory  is  specified,  but  the  velocities  are  so  large 
that  sharp  motions  are  required  near  the  top  of  the  cable,  to  overcome  drag  in  the  lower 
portion.  When  the  frequency  content  of  these  motions  exceeds  that  of  the  retained  Fourier 
series,  the  iterative  procedure  may  be  imstable.  Thus,  the  user  should  work  with  a  model 
of  sufficiently  high  order,  and  keep  only  the  frequency  content  required  for  the  application. 

There  are  a  few  inescapable  numerical  flaws  worth  noting.  The  first  is  the  way  in  which 
i{Xk)  is  computed;  the  operator  i  takes  an  in  the  frequency  domain,  into  a  force, 
also  in  the  frequency  domain.  To  do  the  calculation,  it  is  almost  always  necessary  to  take 
the  inverse  transform  of  Xk,  get  the  time-domain  force,  and  then  transform  back  to  the 
frequency  domain.  The  fault  of  this  procedure  is  that  the  two  transforms  are  generally 
approximations  if  any  of  the  series  is  missing,  and  spurious  components  in  #(X)  could  be 
generated  and  amplified.  Another  point  is  that  the  inversion  of  can  never  be  exact, 

which  leads  to  roundoff  errors  in  13.  Finally,  Gibb's  phenomenon  almost  always  leads  to 
poor  behavior  when  sharp  transients  appear  in  any  component  of  z.  This  response  may  be 
somewhat  alleviated  by  filtering  or  windowing  techniques. 

The  results  of  this  section  required  no  particular  attention  to  the  stability  issue,  and  we 
give  a  full  description  of  the  parameters  for  each  run  to  indicate  the  normal  behavior  of  the 
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algorithm. 


Drill  String  with  Quadratic  Drag 

The  first  example  is  a  simulation  in  which  a  full-scale  drill  string  system^  is  subjected  to 
quadratic  drag,  during  an  arbitrary  in-plane  towfish  trajectory.  Figure  3-13  shows  this  run 
as  a  time  history  of  the  endpoints,  and  gives  a  stroboscopic  view  of  the  cable  shape  evolution 
through  the  first  80  minutes  of  the  run.  The  series  truncation  corresponds  to  a  bandwidth  of 
.0037  radians /second,  and  n  + 1  is  thirteen.  The  algorithm  took  nine  iterations  to  converge 
to  within  an  error  norm  of  10~^. 

The  acausal  nature  of  the  algorithm  is  evident  from  the  figure,  as  the  surface  vessel 
motion  precedes  the  towfish  motion  in  both  of  the  startups  at  2000  seconds  and  7000 
seconds.  We  also  see  Gibb’s  phenomena  clearly  in  the  vessel  trajectory  before  the  large 
reverse  move.  Finally,  the  algorithm  has  successfully  recovered  the  overshoot  concept  of 
Section  3.2;  the  vessel  trajectory  beginning  at  7000  seconds  has  a  300-meter  overshoot  to 
move  the  towfish  back  to  the  origin. 

Inclusion  of  Small  Bending  Stiffness  Terms 

Clearly,  when  the  cable  or  riser  has  significant  bending  stiffness,  the  tridiagonality  of  Equa¬ 
tion  3.2  is  lost.  As  such,  convergence  of  the  iterative  scheme  we  have  given  cannot  be 
guaranteed.  Nonetheless,  it  seems  reasonable  that  one  might  try  to  carry  out  the  inversion 
with  bending  stiffness  anyway,  since  it  is  small  for  many  deep  ocean  systems.  In  fact,  this 
turns  out  to  work  quite  well,  as  is  shown  in  Figure  3-14.  The  parameters  for  the  run  are 
exactly  the  same  as  in  the  previous  case,  except  that  we  have  included  a  bending  stiffness 
which  is  ten  times  the  actual  stiffness  for  the  pipe.  The  bending  moments  were  computed 
using  the  standard  finite- difference  formvilas  for  fourth  derivatives.  Ten  iterations  were 
required  for  this  run. 

The  major  feature  here  is  that  the  stiffness  of  the  pipe  enables  the  algorithm  to  use  the 
weight  more  effectively,  requiring  less  overshoot.  Gibb’s  effect  is  still  visible,  but,  as  before, 
does  not  cause  numerical  difficulties  or  significant  deviation  of  the  towfish. 

^see  Chapter  5  for  its  specifications 
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horizontal  excursion,  meters 


Figure  3-13:  Time  series  (top)  and  60-second  strobed  side  view  (bottom)  of  a  flexible  driU 
string  response  during  preshaping  run. 


horizontal  excursion,  meters 


Figure  3-14:  Time  series  (top)  and  GO-second  strobed  (bottom)  side  view  of  a  stiff  drill 
string  response  dtiring  presbaping  rim. 


Coupled  In-Plane  and  Out-of-Plane  Towing 

As  was  shown  in  the  derivations  of  Section  2.5,  the  zeroth-order  three-dimensional  dynamic 
equations  for  a  cable  are  geometrically  uncoupled.  In  the  ocean,  this  means  that  the  drag 
term  alone  links  in-plane  motions  to  out-of-plane  motions.  The  presence  of  a  coupling 
nonlinearity  between  distinct  tridiagonal  physical  systems  can  be  handled  quite  easily  with 
our  technique.  Since  each  iteration  hxes  the  trajectory  of  one  more  physical  point  on  the 
cable,  we  can  solve  the  two  parallel  problems  separately,  point  by  point.  The  key  is  that 
the  nonlinearity  is  a  simple  function  of  the  local  in-plane  and  out-of-plane  velocities,  which 
are  known. 

Figure  3-15  shows  how  this  procedure  may  be  used  to  compute  a  vessel  trajectory  that 
leads  to  a  preplanned  response  of  the  lower  endpoint  of  the  cable.  For  this  example,  a 
nominal  ARGO/ JASON  system  with  3000  meters  of  cable®  is  to  execute  a  200-meter  radius 
towfish  circle  at  one-half  knot.  The  surface  vessel  trajectory  is  essentially  a  larger  radius 
circle,  with  starting  and  stopping  transients.  The  fact  that  the  ship  trajectory  is  not  round 
is  an  artifact  of  the  series  truncation,  as  before.  Sixteen  components,  for  a  bandwidth  of 
.0049  radians /second,  were  retained  for  the  run,  which  involved  ten  iterations. 

It  can  be  verified  that  the  vessel  trajectories  required  for  rectangular  grid  patterns  come 
out  to  be  essentially  the  same  as  those  used  in  Section  3.2.  As  a  fin2d  note,  in  the  deep 
ocean,  there  are  often  known  currents  which  have  a  significant  effect  on  the  motions  of  the 
cable.  These  currents  can  be  incorporated  into  the  scheme  easily  by  explicitly  including 
them  in  the  time- domain  calculation  of  drag  forces. 

The  Heat  Equation  with  Radiative  Cooling 

Other  problems  can  be  solved  with  the  inversion  method.  This  example  considers  a  situation 
in  which  the  far  end  of  a  long  metal  bar  is  to  be  heated  by  applying  heat  to  the  near  end.  The 
application  might  be  annealing  in  a  foimdry,  where  very  precise  temperature  trajectories 
are  required,  or  a  long  rod  in  space.  We  suppose  that  the  bar  is  long  and  thin  enough  that 
the  one- dimensional  heat  equation  is  valid,  and  we  assume  linear  convective  cooling  of  the 
rod,  as  well  as  radiative  cooling.  We  have,  in  nondimensional  form  [59], 

*see  Chapter  5  for  its  specifications 
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y  meters 


-  T.m««a{*,t))  -  c^T\x,t),  (3.23) 

where  T  is  the  nondimensional  temperature  in  the  rod,  t  is  nondimensional  time,  and  Qc  and 
Or  are  positive  constants  controlling  the  convective  and  radiative  heat  losses,  respectively. 
The  boundary  condition  is 


^  (3.24) 

dt  dt  '  '  ^ 

where  Q  is  the  heat  going  into  the  bar  at  x  =  0.  The  form  of  the  equation  is  familiar;  it 

clearly  can  be  cast  in  a  tridiagonal  finite-difference  approximation,  and  the  nonlinearity  is 

local,  as  required.  An  example  of  the  inversion  method  applied  to  this  problem  is  given  in 

Figure  3-16.  The  parameter  a,,  was  taken  to  be  unity,  and  Or  was  set  to  .01.  The  ambient 

temperature  was  assigned  a  sinusoidal  time  dependence,  and  zero  on  the  temperature  scale 

T  corresponded  with  two  on  the  absolute  scale,  for  computation  of  the  radiation  term.^ 

The  inversion  here  required  seven  iterations,  and  the  number  of  nodes  in  the  model  was 

thirteen;  sixteen  Fourier  components  were  kept,  out  of  a  possible  256. 

The  inversion  scheme  is  able  to  account  for  the  continual  dissipation  of  heat  due  to 
radiation:  the  near  end  temperatures  are  higher  when  the  far  end  temperature  is  high.  It 
also  takes  care  of  the  time-variations  of  the  ambient  temperature:  fluctuations  in  the  near 
end  temperature  are  out  of  phase  with  the  ambient  temperature. 

3.3.5  Preshaping:  Full-Scale  Experimental  Results 

During  the  same  cniise  that  was  described  in  Section  3.2,  tests  were  conducted  to  verify  the 
dynamic  inversion  procedure  in  a  single  plane.  We  made  two  full-scale  tests:  the  results  are 
shown  in  Figures  3-17  and  3-18.  The  top  plot  of  each  set  gives  the  output  of  the  inversion 
program  with  the  towfish  motions  specified  in  advance.  In  the  first  case,  the  clump  weight 
is  required  to  accelerate  from  rest  and  travel  at  constant  speed  for  thirty  minutes,  then 
decelerate  to  a  slower  speed  for  another  thirty  minutes,  then  reverse  and  follow  a  similar 
profile  back  to  the  original  location.  The  inversion  required  six  iterations  for  a  five-node 
plant  model,  with  twelve  frequency  components  retained.  In  the  second  case,  the  clump 
weight  is  only  to  accelerate  to  a  constant  speed  and  then  stop  after  approximately  thirty 

*In  dimensional  units,  zero  on  our  T  scale  would  correspond  to  273  kelvins. 
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minutes.  This  run  also  took  six  iterations  to  converge.  The  time  scales  chosen  for  these 
two  plots  reflect  the  fact  that  the  experimental  runs  were  initiated  at  different  points  in 
the  “total  trajectory”  specified  by  the  inversion;  that  is,  the  second  test  was  given  a  much 
longer  lead-in  time. 

The  experimental  results  are  shown  in  the  lower  portions  of  the  figures.  In  the  first  test, 
only  one  preliminary  oscillation  in  the  ship  motion  was  included,  and  the  clump  weight 
moved  a  significant  distance  in  response.  A  very  dean  movement  from  twenty  to  forty-two 
minutes  elapsed  followed,  and  the  change  in  vdodty  at  forty-two  minutes  was  sharp.  The 
lower-vdodty  forward  motion  was  very  smooth  also.  We  lost  GPS  tracking  at  seventy-two 
minutes  and  were  unable  to  regain  it  until  one-hundred  minutes  time  elapsed.  As  a  result, 
the  latter  part  of  the  nm  was  made  with  a  delay  of  approximately  twenty  minutes.  The 
cable  evidently  assumed  the  proper  shape  agsun  in  time  to  achieve  a  very  sharp  stop  at  142 
minutes  elapsed,  within  seven  meters  of  the  starting  point. 

For  the  second  test,  the  lead-in  oscillations  were  included  in  the  vessel  trajectory.  As 
shown  in  Figure  3- .18,  this  extra  motion  resulted  in  some  movement  of  the  clump  weight, 
but  qualitatively,  the  acceleration  beginning  at  thirty  minutes  elapsed  was  better  than  in 
the  first  run.  The  constant-speed  period  was  quite  smooth,  and  the  deceleration  at  the  end 
of  the  run  was  sharp.  The  final  position  error  of  approximately  forty  meters  can  be  largely 
attributed  to  the  response  of  the  clump  weight  to  oscillations  early  in  the  run. 

A  side  view  of  the  motions  of  the  towfish  for  the  second  preshaping  nm  is  given  in 
Figure  3-19.  In  comparison  with  Figure  3-6  for  the  pure  overshoot  case,  this  nm  is  more 
aggressive,  in  the  sense  that  the  vertical  excursion  of  the  towfish  was  twice  as  large.  This 
difference  is  largely  due  to  the  fact  that  vessel  velocities  for  the  present  nm  were  near  0.7 
meters  per  second,  compared  to  0.5  meters  per  second  in  the  overshoot  case.  Nonetheless, 
the  vertical  excursions  should  not  be  a  problem  for  most  winches  to  keep  up  with,  as  the 
maximum  rate  of  altitude  change  for  this  nm  is  under  0.25  meters  per  second. 


3.4  Summary 

The  results  of  this  chapter  indicate  that  the  idea  of  input  preshaping  is  an  effective  way 
both  to  improve  the  speed  of  response  in  simple  point-to-point  moves  of  the  towfish,  and  to 
execute  complex,  coupled  trajectories.  The  full-scale  experimental  data  show  that  effective 
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Figure  3-18:  Analytical  trajectory  (top)  and  experimental  results  (bottom)  for  the  second 
full-scale  test. 
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Figure  3*19:  Measured  heave  response  during  the  second  full-scale  preshaping  test. 


dynamic  positioning  of  the  surface  vessel  plays  a  crucial  part;  the  runs  could  have  been 
cleaner  if  the  DP  system  had  the  capability  to  follow  predefined  trajectories,  for  example. 
At  the  same  time,  the  inherent  safety  of  this  technique  has  been  illustrated  through  the 
navigation  failure  in  the  first  preshaping  test.  Since  the  procedures  involve  no  feedback,  such 
a  failure  means  only  that  the  plant  will  deviate  from  the  desired  track,  and  return  smoothly 
to  it  after  navigation  is  restored.  Finally,  the  vertical  excursions  of  the  towfish  during 
these  tests  were  slow  enough  that  typical  winches  should  have  no  difficulty  in  maintaining 
constant  altitude. 
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Chapter  4 


Tools  for  the  Cable  Feedback 
Control  Problem 


4.1  Introduction 

This  chapter  considers  feedback  control  for  positioning  of  a  towed  cable  system.  The  prob¬ 
lem  is  important  because  many  deep-ocean  operations  with  towed  vehicles  require  accurate 
positioning  in  the  presence  of  unknown  disturbances  such  as  currents.  Further,  if  the  depth 
is  large,  and  the  drag  forces  are  large,  then  the  plant  is  clearly  infinite-dimensional  and 
nonlinear.  These  properties  make  the  problem  difficult,  as  the  guarantees  built  into  many 
control  techniques  cannot  be  valid.  In  Chapter  2,  an  attempt  was  made  to  justify  each 
step  of  the  model  simplification  process.  Now,  we  investigate  performance  and  robustness 
guarantees  with  respect  to  the  reduced-order  models.  We  build  upon  a  great  deal  of  work 
by  other  researchers,  and,  as  such,  a  large  part  of  this  chapter  covers  known  material. 

4.2  Overview 

Since  the  internal  state  of  a  cable  (between  the  endpoints)  clearly  has  a  dominant  role  in 
the  dynamic  response  at  the  endpoint,  full-state  feedback  schemes  require  adequate  sensors 
or  state  estimation.  There  are  at  least  several  conceivable  ways  in  which  state  measurement 
could  be  carried  out  for  a  deep-water  cable  system.  First,  additional  transceivers  could  be 
affixed  to  the  cable;  given  the  calculations  of  Section  2.6.2,  perhaps  only  a  few  would  be 
reqtiired.  Alternatively,  measurement  of  the  cable  angle  at  the  surface  could  be  used  in  an 
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energy-based  feedback  law  (e.g.,  [92]).  The  idea  would  be  to  move  the  vessel  in  such  a  way 
that  energy  is  always  being  extracted  from  the  cable.  Both  of  these  approaches,  however, 
are  difficult  to  envision  in  heavy  weather — the  cable  angle  would  be  an  extremely  noisy 
measurement,  and  the  attachment  of  transceivers  on  the  cable  would  be  very  dangerous. 
For  these  reasons,  we  assume  in  this  chapter  that  the  only  available  navigation  is  of  the 
towfish  and  of  the  siirface  vessel.  With  typical  long-baseline  acoustic  navigation  systems 
and  GPS,  the  noise  variance  of  both  measurements  is  approximately  five  meters  or  less. 

The  use  of  estimators  in  closed-loop  control  has  been  a  central  issue  in  the  literature 
for  at  least  two  decades.  For  the  linear  case,  the  counterexample  of  Doyle  [36]  showed 
that  arbitrary  coupling  of  Kalman  filters  with  linear  quadratic  regulators  (LQR)  carries  no 
firm  robustness  guarantees,  even  though  the  robustness  properties  of  each  part  are  good; 
such  guarantees  had  been  h3rpothesized  by  Safonov  and  Athans  [84].  The  existence  of  this 
counterexample  was  sufficient  to  motivate  the  Loop  Transfer  Recovery  (LTR)  procedure, 
which  has  been  a  widely  favored  version  of  the  LQG  approach  for  a  decade. 

To  date,  there  is  no  accepted  analogue  for  nonlinear  systems,  although  Safonov  gave  a 
nonlinear  separation  principle,  which  we  state  here  (see  Appendix  A  for  definitions): 

Theorem[82]:  Suppose  a  nonlinear  feedback  system  is  closed-loop  bounded  (finite-gain 
stable)  with  the  feedback  u  =  — p(i),  where  <7  is  a  causal  nonlinear  d}mamical  operator  with 
finite  incremental  gain.  Let  £  be  an  estimate  of  x  which  is  nondivergent  (with  finite  gain). 
Then  the  feedback  system  with  x  replacing  x  in  the  control  law,  is  closed-loop  bounded 
(finite-gain  stable). 

The  idea  is  that  if  a  nondivergent  estimator  can  be  obt2Lined  for  the  system,  then  one  can 
implement  a  possibly  nonlinear  control  law  based  on  exact  state  measurement  and  expect 
the  closed-loop  system  to  perform  well.  As  with  linear  separation,  however,  there  is  no  ac¬ 
commodation  for  uncertainty  in  the  plant  model.  To  get  around  this  shortcoming,  Grunberg 
and  Athans  [50], [51]  applied  the  idea  of  loop  transfer  recovery  to  nonlinear  systems;  their 
methodology  was  termed  Loop  Operator  Recovery  (LOR).  Like  LTR,  the  LOR  has  dual 
formulations,  one  of  which  recovers  the  controller  loop  and  the  other  the  estimator  loop. 
Unfortunately,  the  more  useful  estimator  loop  recovery  procedure  has  been  proven  only  for 
systems  in  controUer  and  observer  form  (e.g.,  robots).  Recovery  of  tbe  controDer  loop  is 
much  easier  to  achieve  theoretically,  although  it  puts  severe  limitations  on  the  performance 
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and  robustness  specifications. 

The  control  techniques  that  are  considered  in  this  chapter  include  nonlinear  model-based 
observation,  and  the  design  of  nonlinear  regulators,  for  direct  use  with  Grunberg’s  control 
loop  recovery  procedure.  Thus,  we  first  outline  what  L0&  entuls,  and  detail  some  prop¬ 
erties  of  various  estimators  (nominal  nondivergence,  parameterization)  and  control  loops 
(robustness  and  performance).  Then  the  primary  contribution  of  the  chapter  is  developed: 
a  pr&ctical,  approximately  optimal  nonlinear  feedback  law  for  high-order  plant  models. 

4.3  Loop-Operator  Recovery 

In  this  section,  we  describe  more  specifically  the  LOR  technique  that  is  to  synthesize  nondi- 
vergent  estimators  with  the  LQR  or  the  nonlinear  optimal  controller.  Along  with  listing  the 
relevant  theorem,  the  practical  shortcomings  of  the  approach  are  discussed,  to  be  revisited 
experimentally  in  Chapter  5. 

4.3.1  The  LQG/LTR  Technique 

What  follows  is  a  review  of  loop  transfer  recovery  [39],  [94],  given  because  such  a  discussion 
is  helpful  to  the  uninitiated,  and  because  the  nonlinear  LOR  technique  of  the  next  section 
is  in  complete  analogy  with  LTR.  Our  discussion  is  brief  because  one  need  not  use  LTR 
to  solve  the  linear  SISO  problem  associated  with  our  cable  system,  and  because  LTR  has 
been  the  topic  of  countless  papers  and  books  in  the  last  decade.  LTR  takes  two  forms;  one 
hinges  upon  the  good  properties  of  the  LQR  loop,  and  the  other  takes  advantage  of  the 
similar  properties  of  the  Kalman  filter  loop.  As  noted  previously,  we  will  be  recovering  the 
control  loop  only. 

The  first  step  is  to  carry  out  the  LQR  design  procedure  for  the  plant  {si  — 
with  Q  =  C^C  and  p  chosen  to  put  the  crossover  at  an  appropriate  frequency.^  This 
gives  a  feedback  gain  matrix  G,  and  the  control  loop  G{sl  -  A)~^B  has  a  good  shape 
for  robustness  against  modeling  errors  and  disturbances  at  the  plant  input.  Although  the 
command-following  and  output  disturbance  rejection  properties  of  such  a  loop  cannot  be 
manipulated  explicitly,  these  designs  are  usually  acceptable  nonetheless,  at  least  for  for 
minimum-phase  plants  [4].  The  control  loop  transfer  function  ideally  is  very  large  at  low 

‘The  notation  is  established  in  Sections  4.4  and  4.5. 
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frequencies,  passes  smoothly  through  the  crossover  region,  and  attenuates  rapidly  at  high 
frequencies. 

To  recover  tL  .a  target  loop  while  using  a  state  estimator,  we  choose  L  =  B,  and  let  p  be 
a  free  parameter  in  the  filter  design.  It  can  be  shown  that  liin^_o  =  BW,  where  W  is 
an  orthononnal  matrix,  and  eventually  that  as  /i  — »  0,  K{s)P(a)  — »  G{sl  -  A)~^B.  Here, 
H  is  the  filter  gain,  K{s)  is  the  control  transfer  matrix,  and  P(s)  =  C{sl  -  A)~^B  is  the 
plant.  Thus,  we  recover  all  of  the  robustness  and  performance  properties  that  accompany 
the  LQR  loop. 

4.3.2  Nonlinear  Extension — Loop  Operator  Recovery 

The  idea  of  nonlinear  loop-operator-recovery  (LOR)  is  completely  analogous  to  the  LTR 
results  above.  We  present  here  only  one  of  three  variations,  since  the  other  two  (LOR  at 
the  plant  output  and  formal  loop  shaping)  require  conditions  that  are  overly  restrictive  [50]. 

The  method  calls  for  recovery  of  the  controller  loop.  Thus,  the  control  law  itself  can 
be  designed  with  any  technique,  but  the  main  requirement  is  that  it  yield  a  loop  with 
good  robustness  and  performance  properties.  If  a  nonlinear  optimal  solution  is  used,  these 
properties  have  been  found  to  be  very  good  (see  Section  4.5.2).  As  in  the  LTR  situation, 
we  then  use  a  parameterized  filter  such  that  as  the  sensor  noise  intensity  goes  to  zero,  the 
full  closed-loop  system  with  the  observer  emulates  the  pure  feedback  system  with  perfect 
measurements.  The  crucial  theorem  is  as  follows: 

Theorem[51]:  Let  the  filter  gain  be  parameterized  by  p  such  that  lim^_o  =  BW, 

where  W  is  some  orthonormal  matrix.  If  H  is  linear  then  lim^_o(-'^^)(~'P)  =  GiB,  where 
P  is  the  nonlinear  pliuit  operator,  K  is  the  control  operator,  and  the  nonlinear  plant  is  iB. 


In  this  theorem,  the  control  law  G  may  be  nonlinear,  time- varying,  or  could  be  a  dynamical 
system  itself.  The  generality  of  this  result  is  very  attractive,  as  it  allows  any  nominally 
stabilizing  controller  to  be  used,  and  admits  a  wide  range  of  observers.  In  aU  cases,  nominal 
nondivergence  of  the  estimator,  and  stability  of  the  control  law  must  be  maintained. 
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4.3.3  Recovering  the  Regulator  Loop 

It  is  important  to  reiterate  that  this  technique  recovers  not  the  filter  loop,  but  the  control 
loop.  Thus,  in  the  same  way  that  LQR-loop  recovery  in  the  linear  case  is  restrictive,  the 
situation  here  is  not  ideal.  Figiure  4-1  gives  a  simple  explanation  of  why  recovering  the 
control  loop  is  so  restrictive.  Part  a)  gives  the  standard  regulator  loop,  which,  with  LQR 
or  nonlinear  feedback  designs,  has  good  properties  broken  at  the  plant  input  (x).  Part 
b)  is  an  expanded  version  of  a);  we  have  broken  the  state  vector  into  observed  (y)  and 
unobserved  (z,)  parts.  With  the  loop-breaking  point  (x)  at  the  same  spot  in  the  overall 
loop,  we  see  that  the  properties  of  the  loop  broken  at  the  point  (xx)  are  not  the  same. 
The  latter  location  is  important  in  terms  of  output  disturbance  rejection  and  command 
following.  In  short,  we  cannot  get  direct  access  to  the  output  y,  the  variable  upon  which 
output  disturbances,  sensor  noise,  and  reference  signals  act  directly.  Also,  we  must  place 
our  modeling  uncertainty  at  the  plant  input,  along  with  the  disturbances. 

In  the  absence  of  a  better  approach,  we  have  to  settle  for  designing  the  loop  to  meet 
robustness  conditions  at  the  plant  input,  and  for  taking  whatever  command-following  and 
output  disturbance  rejection  characteristics  come  with  the  loop.  For  minimum-phase  linear 
systems,  it  has  been  found  that  the  performance  can  be  still  quite  good  [4],  and  we  expect 
a  similar  trend  for  the  nonlinear  case. 
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4.4  State  Estimation 


4.4.1  Kalman  Filter 

We  outline  the  well-known  linear  observation  problem  to  start,  since  most  nonlinear  ob¬ 
servers  are  similar  in  form.  The  general  n-state  linear  system  is  defined  as 

=  Ax  +  Bu  -I-  Li  (4.1) 

=  Cx  6^ 

where  x  is  the  n  x  1  state  vector,  u  is  the  control  input,  and  y  is  the  observed  output.  A 
is  the  system  matrix,  B  the  input  gain  matrix,  C  the  output  gain  matrix.  The  system  is 
affected  by  process  noise  i,  which  acts  through  the  gain  matrix  L  and  has  intensity  E,  and 
also  by  meastirement  noise  0  of  intensity  G. 

The  Kalman  Filter  (KF)  takes  the  form  [44] 

^  =  Ax  +  Bu  +  E(y-Cx),  (4.2) 

at 

for  which  the  filcer  gain  matrix  H  is  computed  off-line  as  H  —  SC^G~^ .  The  error  covariance 
matrix  £  is  the  symmetric  positive  definite  solution  of  the  filter  algebraic  Riccati  equation: 

AL  -f-  -h  iEX^  -  E<7^0-^CS  =  0.  (4.3) 

The  solution  £  is  guaranteed  to  exist  if  [A,X]  is  stabilizable  and  [A,C]  is  detectable;  the 
standard  technique  today  for  solving  algebriuc  Riccati  equations  involves  a  Schiir  vector 
decomposition  [61]. 

The  Kalman  filter  is  known  to  be  nominally  nondivergent,  and  to  be  robust  against 
uncertainties  at  the  plant  output  [39].  However,  it  does  not  provide  nondivergent  estimation 
for  nonlinear  plants. 

4.4.2  Extended  Kalman  Filter 

Now  we  consider  t:te  generic  nonlinear  system 


dx 

dt 

y 
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dx 

di 

y 


f{x)  +  +  If 

Cz  +  6, 


(4.4) 


where  /(•)  is  a  nonlinear  operator,  and  the  other  variables  are  as  described  for  Equation  4.1. 

Solving  the  iiill  optimal  estimation  problem  [30]  for  nonlinear  systems  involves  the  real¬ 
time  propagation  of  a  partial  differential  equation,  and  the  memory  and  computational 
requirements  of  such  a  technique  are  usually  prohibitive,  althougih  the  estimate  provided  is 
indeed  nondivergent.  These  practical  difficulties  lead  to  suboptimal  approximations,  one  of 
which  is  the  Extended  Kalman  Filter  (EKF).  The  EKF  was  developed  as  a  direct  extension 
to  the  Kalman  Filter,  and  uses  a  nonlinear  model  with  a  time-varying  filter  gain  matrix 
[44]: 


~  =  m^Bu-^H(y-Cz), 

H  =  EC^0-\ 

~  =  V/(i)E -I- EV/(i)^ XEI^  -  EC^0-^CS. 
at 


(4.5) 

(4.6) 

(4.7) 


The  parallels  between  these  equations  and  those  of  the  Kalman  filter  are  obvious.  For  a 
long  time,  there  was  no  theoretical  justification  for  the  EKF,  although  users  recognized 
heuristicaUy  that  it  had  very  good  properties.  Grunberg  [50]  then  showed  that  the  EKF  is 
in  fact  nondivergent  under  very  broad  conditions.  In  particular,  we  cite  the  following  result: 

Theorem  [51]:  Suppose  |  V/(r)  \<  B/  <  oo  *,  and  that  Ea  is  positive  definite.^  Under 
these  conditions,  if  there  exists  any  nondivergent  estimator  for  the  nonlinear  system,  then 
the  EKF  is  nondivergent. 


Given  that  the  EKF  provides  nondivergent  estimates  for  a  wide  range  of  nominal  nonlin¬ 
ear  plants,  one  can  investigate  its  other  properties  so  that  LOR  can  be  carried  out.  For  the 
current  work,  the  asymptotic  behavior  of  the  gain  matrix  as  the  sensor  noise  goes  to  zero  is 


^This  may  be  relaxed,  leading  to  small-signal  results. 
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most  important,  since  we  will  be  recovering  the  regulator  loop,  to  be  defined  in  Section  4.5. 
With  0  =  /il,  the  condition  we  need  is  as  follows: 

^  =  XW(0,  (4.8) 

where  W  is  an  orthonormal,  time-varying  matrix.  The  proof  that  this  is  the  case  for  the 
KF  and  the  CGEKF  of  the  next  section  is  wdl-known  [60],  but  the  EKF  presents  more  of  a 
difficulty,  since  the  propagation  of  £  is  obviously  dependent  on  the  state  estimate.  To  date, 
the  asymptotic  properties  of  the  EKF  are  known  only  locally,  that  is,  for  initial  £  within  a 
small  neighborhood  of  the  origin  [106].  Nonetheless,  it  may  be  numerically  verified  that  for 
many  Lipschitz  nonlinearities,  the  asymptotic  behavior  we  need  holds  nonetheless  for  the 
EKF. 

The  EKF  loop  does  carry  some  robustness  properties  also;  they  are  analogous  to  those 
of  the  Kalman  filter  loop,  and  have  been  discussed  by  Grunberg  [50].  The  cost  for  all  of 
the  good  attributes  of  the  EKF  is  that  it  is  somewhat  laborious  to  propagate  because  the 
covariance  matrix  £(t)  has  to  be  integrated  forward  in  time.  Although  modem  computers 
make  it  easier  to  use  the  EKF  in  real-time,  there  are  many  situations  in  which  these  demands 
are  simply  too  great. 

4.4.3  Constant-Gain  Extended  Kalman  Filter 

The  Constant-Gain  Extended  Kalman  Filter  (CGEKF)  [84],  is  a  simplification  of  the  EKF, 
as  it  avoids  costly  time  integration  of  the  covariance  matrix  £.  In  short,  the  CGEKF  uses 
the  known  nonlinearities  in  propagating  z,  but  with  a  constant  filter  gain  matrix  H.  H  is 
usually  chosen  based  on  a  nominal  value  of  V /(z),  which  we  call  A.  In  a  sense,  the  CGEKF 
is  a  hybrid  between  the  EKF  and  the  KF,  as  the  filter  g^  is  computed  once  according  to 
Equation  4.3,  and  the  state  estimate  propagates  as  in  Equation  4.5. 

A  condition  for  the  nondivergence  of  the  CGEKF  is  the  following: 

Theorem[84]:  If  uniformly 

(A  -  V/(z))£  -f-  +  £C^0-^C£)  >  0,Vz3i2"  (4.9) 

A 

then  the  CGEKF  is  nondivergent  with  finite  giun. 
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A  similar  condition  has  been  provided  by  Tam  and  Rasis  [95] .  Needless  to  say,  these 
conditions  must  be  checked  for  a  range  of  operating  points,  since  the  robustness  measures 
are  with  respect  to  the  nominal  linf'Ar  design  plant.  If  an  uncertain  linear  element  is  placed 
at  the  plant  output,  the  condition  listed  above  yields  iniinite  upward  gain  margin,  50  percent 
gain  reduction  tolerance,  and  at  least  60  degrees  of  phase  margin  in  every  output  channel. 
Thus,  the  CGEKF  can  be  used  in  loop-operator  recovery  of  the  controller  loop,  but  the 
nondivergence  must  be  checked  numerically  throughout  the  operating  range. 

A  critical  shortcoming  of  the  CGEKF  is  that  there  is  no  robustness  to  variations  in  the 
nonlinear  plant.  The  model  in  Equation  4.5  must  be  exactly  right.  As  a  result,  the  dual 
LOR  approach  is  pointless  with  the  CGEKF. 

4.4.4  Other  Nonlinear  Observers 

This  list  of  approaches  (from  Section  4.4.1  on)  is  by  no  means  exhaustive,  for  there  exist 
many  other  nonlinear  observation  techniques.  The  ones  above  are  listed  because  they  are 
among  the  simplest,  and  because  they  have  some  guaranteed  nondivergence  properties. 
This  nondivergence  with  respect  to  a  nominal  nonlinear  plant  is  crucial  in  the  application 
of  LOR.  For  some  other  observer  ideas,  the  reader  is  referred  to  the  survey  paper  by  Misawa 
and  Hedrick  [68],  and  the  references  therein,  as  well  as  papers  by  Thau  [96],  Slotine  et.  al 
[91],  and  Bestle  and  Zeitz  [17]. 

4.5  Full-State  Feedback 

We  now  consider  the  complementary  problem  to  observation:  feedback  assuming  perfect 
measurements.  As  in  the  LTR  technique,  the  eventual  properties  of  the  nonliaear  closed 
loop  depend  on  either  the  estimator  or  the  control  law.  Since  our  strategy  is  to  recover  the 
control  loop,  we  pay  special  attention  to  the  loop  properties  associated  with  these  control 
laws. 


4.5.1  Nonlinear  State  Feedback 

The  imperturbed  nonlinear  system  is  defined  as  follows: 

^  =  /(.)  +  Bu.  (4.10) 
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Our  problem  is  to  find  the  control  u  =  — p(z)  such  that  the  cost  function  J  is  minimized: 


J  =  ^  +  u^Ru]dt.  (4.11) 

2  Jo 

Of  course,  different  cost  functionals  may  be  used,  involving  more  complex  combinations  of 
X  and  u,  and  terminal  constraints.  We  will  keep  the  quadratic  form,  however,  because  it  is 
a  requisite  for  the  main  theorem  we  wish  to  use.  Another  constraint  is  that  4.10  is  linear 
in  the  input;  this  point  also  is  important  for  our  method,  and  is  consistent  with  the  models 
in  Section  2.7. 

There  are  two  parallel  approaches  for  solving  this  optimization  problem;  dynamic  pro¬ 
gramming  and  Pontryagin’s  Maximum  Principle. 


Dynamic  Programming 

The  dynamic  programming  approach  is  due  to  Bellman  [14];  in  it,  one  defines  the  optimal 
return  function  as  the  cost  associated  with  the  optimal  trajectory,  beginning  at  some  inter¬ 
mediate  time  t.  For  example,  with  the  cost  functional  given  by  Equation  4.11,  the  optimal 
return  is 

1 

V{x,t)  =  min  -  I  [x^Qx  RvL]dt.  (4.12) 

“  2  Jt 

The  Hamilton- Jacobi-Belbnan  (HJB)  equation  provides  the  mechanism  for  minimizing 
V(x,t)  at  each  point  along  the  way: 


dV(x,t) 

dt 


-f-  min 


1  Tq  I  ,  dV{x,t) 

-*  RU+-J— 


(/(r)-hBu)]  =0.  (4.13) 


with  the  terminal  condition  that  V{x{T),T)  =  0.  The  minimization  is  made  by  setting  the 
derivative  of  the  bracketed  expression  with  respect  to  u,  equal  to  zero.  This  step  gives  the 
control  u  as  a  ftmction  of  V{x,  t),  namely 


u  = 


dVjx^tf 

dx 


which  leads  to  another  form  of  the  HJB  equation: 
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(4.14) 


dV{x,t) 

dt 


ldV{x,t)„„_,  rdV(x,tf 

^  '"IT' 


dv{x,t) 

dx 


f{x)  =  0. 


Inspection  shows  that  we  have  a  partial  differential  equation  to  solve,  and  since  V(x(T),  T)  = 
0,  it  has  to  be  solved  backwards  in  time.  In  the  limit  as  T  — »  oo,  however,  the  partial  time 
derivative  of  V(z,t)  can  be  set  to  zero,  as  we  expect  it  to  approach  some  constant  value 
toward  t  =  0.  The  remaining  terms  in  this  case  form  an  algebraic  matrix  eqiiation  which 
has  to  be  solved  for  at  all  points  on  the  optimal  trajectory  x,.  In  practice,  one  could 

set  up  a  grid  in  some  region  of  the  state  space,  and  solve  Equation  4.14  a  priori  for 
at  each  point,  storing  the  results.  The  construction  of  such  a  “look-up”  table  would  require 
a  way  to  store  and  access  these  data  effectively  [40]. 

Alternatively,  Lukes  [65]  and  Willemstein  [103]  developed  a  technique  for  solving  the 
HJB  equation  by  a  power  series  expansion  about  the  origin.  Given  analycity  of  /,  they 
expanded  the  control  u  and  the  optimal  retiim  function  T^(z,  t),  and  with  the  HJB  equation 
and  a  similar  relation,  showed  that  local  convergence  of  the  control  law  holds.  This  technique 
has  been  used  effectively  by  many  authors,  including  Dabbous  and  Ahmed  [32]  and  Dwyer 
and  Sena  [41].  As  with  all  perturbation  techniques,  sometimes  an  anal3rtical  solution  can 
be  achieved,  but  the  numerical  advantages  are  also  significant  in  their  own  right. 


Minimum  Principle 

Another  way  to  formulate  the  optimal  control  problem  involves  Pontryagin’s  Minimum 
Principle.  It  begins  with  the  definition  of  the  Hamiltonian: 

H{x,  u.  A,  t)  := 

where  A  is  the  n  x  1  costate  vector.  The  costate  propagates  according  to 

dX  dJI{x,u,X,t)^ 
di~  ^ 


dX 

dt 


dx 


A-Qz, 


(4.16) 


with  the  terminal  condition  that  A(T)  =  0.  The  minimtun  of  the  Hamiltonian  is  found  by 
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differentiating  with  respect  to  u,  giving 


du 

u  =  (4.17) 

Thus,  the  link  between  the  dynamic  programming  approach  and  the  TniniwiiiTTi  principle 
approach  is  that  X  =  ,  where  A  is  to  be  propagated  backwards  in  time  like  V. 

In  the  same  spirit  of  the  expansions  of  Lukes  and  Willemstein,  Yoshida  and  Loparo  [105] 
devised  a  perturbation  technique  to  solve  the  optimal  control  problem  in  the  context  of  the 
iniTiiTtiuTn  principle.  However,  these  authors  were  able  to  show  not  just  local  convergence 
of  the  power  series  for  the  control  u,  but  global  convergence.  Starting  from  this  expansion, 
we  will  describe  a  general  framework  for  implementing  approximately  optimal  nonlinear 
controllers  in  Section  4.6. 

4.5.2  Properties  of  Optimal  Nonlinear  Control  Loops 

Optimal  nonlinear  control  loops  have  nominal  stability,  tunable  performance,  and  good 
robustness  properties.  First  we  define  the  operator  i  to  take  Bu  to  z  as  4  :=  (s  — 
where  s  is  the  derivative  operator.  With  this  notation,  the  nonlinear  system  of  Equation  4.5 
can  be  written  simply  as  y  =  CiBu. 

To  begin,  the  “return  difference  condition”  of  optimal  controllers  [70],  [51]  is  that 

l|(J  +  G$H)-Ml2.r<l.  (4.18) 

where  the  control  law  k  u  =  —G{x).  Simple  manipulations  from  this  point  show  also  that 

II  +  GiB)-^  ||2,^<  2,  and  (4.19) 

||(J  +  (G«B)-)u||,..  1 

-2  '  ’ 

Equation  4.18  clearly  implies  that  the  system  is  nominally  closed-loop  stable.  In  addition. 

Equation  4.20  provides  stability-robustness  as  in  the  LQR  case;  for  example,  the  presence 
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Figure  4-2:  Two  model  uncertainty  structures:  a)  pre-multiplicative,  and  b)  pre-division. 

of  a  multiplicative  imstnictured  imcertainty  A  (see  Figure  4-2)  at  the  plant  input  can  be 
tolerated  so  long  as  ||  A  ||2,r<  |  [50].  Similarly,  through  Equation  4.18,  a  division  error 
with  II  A  ||2,r<  1  can  be  tolerated.  These  margins  are  identical  to  those  of  LQR  designs: 
50  percent  gain  reduction  tolerwce,  infinite  upward  g^  margin,  and  60  degrees  of  phase 
margin  in  all  input  channels  [99].  However,  it  is  noteworthy  that  these  margins  are  all  to 
be  taken  with  respect  to  the  nominal  nonlinear  pl2mt. 

Good  performance  is  also  provided  by  the  nonlinear  optimal  regulator,  and  as  in  the 
LQR,  the  primary  handle  on  performance  is  the  control  penalty  p.  A  standard  result  is  that 
stability  with  large  gain  implies  good  tracking  performance: 

Theoretn[35]:  Let  R  be  a  set  of  reference  trajectories  r  of  interest.  If 

II  (/  -H  l|2.r«||  r  ||2,r,Vr3R,  then  ||  r  -  u  ||2.r«||  r  Ib.r,  Vr3R. 

In  linear  systems,  the  reference  signals  in  R  are  generally  of  low  frequency.  For  nonlinear 
systems,  it  is  not  so  clear  what  is  in  R,  although  intuitively  the  signals  still  should  be  slow. 
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This  theorem  is  consistent  with  the  fact  that  integral  action  drives  the  steady-state  errors 
to  zero,  even  in  nonlinear  systems  [34].  Furthermore,  the  optimal  nonlinear  control  loop  is 
insensitive  to  disturbances  at  the  plant  input,  although  the  controller  is  obviously  unable 
to  distinguish  between  the  disturbance  and  a  change  in  the  reference  signal. 

The  optimal  controller  also  possesses  a  parameterization  similar  to  the  LQR  cheap 
control  result,  which  makes  it  possible  to  recover  the  nonlinear  filter  loop  in  the  dual  LOR 
procedure  [51].  Again,  this  thesis  considers  only  recovery  of  the  control  loop. 

4.5.3  LQR  Designs 

We  now  consider  the  application  of  a  constant  linear  feedback  law,  as  found  through  the 
LQR  procedure.  Namely,  we  let  the  feedback  matrix  be  67  =  where  K  is  the 

symmetric  positive  definite  solution  to  the  algebraic  Ricatti  equation 

0  =  KA  +  A'^K  +  Q-  -KBB^K. 

P 

Safonov  and  Athans  gave  a  condition  for  closed-loop  finite-gain  stability  with  nonlinear 
plants  [83];  we  must  have  uniformly 

K(A  -  Vf{z))  +  liQ  +  ^KBB'^K)  >  0,  VarBR".  (4.21) 

2  p 

This  matrix  is  generally  computed  with  respect  to  the  nominal  state  linearization,  A  = 
^/(zo)-  The  condition  may  be  massaged  to  reveal  that  the  LQR  control  law  has  the  50 
percent  gain  reduction  tolerance,  infinite  upwards  gain  margin,  and  60  degrees  of  phase 
margin  in  all  input  channels.  In  contrast  to  the  margins  of  optimal  nonlinear  control  laws, 
however,  these  apply  to  the  nominal  linear  plant,  so  we  are  likely  to  lose  this  guarantee  if 
the  plant  ventures  far  away  from  the  nominal  operating  point. 


4.6  Approximate  Nonlinear  Optimal  Feedback  by  Pertur¬ 
bations 

The  above  developments  of  other  researchers  indicate  that  the  optimal  nonlinear  feedback 
law  possesses  attractive  properties  in  terms  of  the  recovered  control  loop.  These  include 
tunable  performance  and  robustness  at  the  plant  input,  with  respect  to  a  nonlinear  nominal 
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plant,  and  so  synthesis  of  such  a  control  law  with  a  nondivergent  observer  such  as  the  EKF 
can  be  expected  to  yield  good  results  through  the  LOR  procedure.  Thiis,  in  the  current 
section  we  present  a  framework  for  the  analytical  approximaticm  of  the  optimal  nonlinear 
feedback  law. 

A  feedback  controller  for  the  cable  problem  of  interest  needs  to  have  several  specific 
attributes  in  order  to  be  tractable.  First,  it  is  apparent  that  a  nondynamical  feedback 
law  would  be  advantageous.  Given  the  computational  effort  that  is  likely  to  be  required 
for  nonlinear  observation,  it  makes  sense  not  to  have  two  dynamical  systems  running  in 
software.  Also,  the  control  law  should  be  computed  in  real-time,  without  resorting  to  large 
lookup  tables.  Solving  the  HJB  equation  a  priori  on  a  grid  may  be  effective  for  low-ordered 
systems,  but  can  be  very  difficiilt  when  the  plant  order  is  high,  as  in  our  cable  models. 

The  perturbation  method  of  Yoshida  and  Loparo  [105]  is  suitable  in  that  it  leads  to 
nondynamical,  real-time  control  laws.  The  basic  idea  is  to  construct  auxiliary  vectors  which 
are  nonlinear  functions  of  the  state  vector,  and  then  find  the  appropriate  matrices  that  take 
these  new  vectors  linearly  into  the  command  signal.  As  will  be  seen  in  the  sequel,  however, 
these  matrices  are  nontrivial  to  find;  in  the  literature,  one  finds  only  examples  of  second- 
and  third-order  systems,  for  which  these  matrices  can  be  worked  out  by  hand.  Clearly,  for 
a  nonlinear  regulator  to  control  a  cable  system,  we  must  achieve  the  generality  to  construct 
the  necessary  matrices  for  plants  of  arbitrary  order.  This  generalization  of  the  pertiirbation 
method  is  presented  in  this  section. 

4.6.1  Foundations  of  the  Perturbation  Approach 
Suppose  f{x)  is  analytic  and  can  be  written  as 


/(*)  =  (4.22) 

k=2 

where  A  is  the  Jacobian  of  /(*)  evaluated  at  the  origin.  The  vectors  xl*"!  contain  terms  of 
order  z^,  and  are  to  be  specified  later,  although  it  is  apparent  that  the  actual  nonlinearity 
should  be  modeled  as  a  polynomial  series.  We  also  assume  that  the  costate  A  is  of  the  form 

A  =  f;PfczW.  (4.23) 

fc=i 

This  formulation  is  completely  dual  to  the  Maximum  Principle  derivation  of  the  LQR,  and 
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gives  rise  to  a  perturbation-type  technique  for  successively  finding  the  matrices  Pk.  Yoshida 
and  Loparo  showed  that  this  particular  expansion  has  the  following  property: 

Theorem  [105]:  Suppose  that  the  xmforced  system  of  Equation  4.5  does  not  have  a  finite 
escape  time.  If,  in  addition,  /  is  an  analytic  function  of  the  state,  /(O)  =  0,  and  the  set 
[i4,  B,  C]  is  controllable  and  observable,  then  the  optimal  regulator  given  by 

u  =  ^  (4.24) 

k=l 

converges  absolutely  and  uniformly  on  [0,  oo],  for  any  z(0)Bi2''.  This  regulator  also  satisfies 
the  HJB  equation. 

In  the  same  paper,  a  necessary  and  sufficient  condition  for  the  optimal  regulator  to  have 
a  finite  series  expansion  was  derived,  although  it  is  usually  difficult  to  verify. 

We  now  outline  the  relevant  expanded  equations,  to  generate  the  optimal  control  law 
above.  With  an  infinite  horizon,  ^  =  0,  for  all  k.  Then  substituting  Equations  4.22 
and  4.23  into  Equation  4.16  gives 


0  =  Qr  -f-  f;  f;  4-  f)  + E  -  jif  f;  p**W).  (4.25) 

fc=2  fc=l  fc=l  fc=2  fc=l 

where  K  :=  BR~^B^.  The  subscript  “r”  denotes  the  first  partial  derivative  with  respect 
to  X,  defined  such  that  /*(*  =  0)  =  A.  On  collecting  terms  of  like  powers,  we  obtain  the 
first-order  equation  as  follows: 

Q  -I-  A^Pi  -I-  PiA  -  PiBR-^B^Pi  =  0.  (4.26) 

The  reader  will  notice  that  this  is  just  the  matrix  algebraic  Riccati  equation.  The  LQR 
with  a  plant  linearized  at  the  origin  can  therefore  be  seen  as  an  estimate  of  the  loop  shape 
before  the  nonlinear  terms  are  added. 

4.6.2  Obtaining  the  Higher-Order  Matrices 

From  this  point  onward,  we  present  the  development  specific  to  a  plant  with  a  symmetric 
cubic  polynomial  approximation  (like  quadratic  drag),  as  a  special  case  of  the  perturbation 
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expansion.  Since  this  means  that  the  even-powered  terms  disappear,  and  that  we  will 
not  consider  terms  higher  than  third-order,  some  of  the  generality  is  lost  at  this  point. 
Nonetheless,  the  derivations  will  indicate  the  nature  of  the  calculations  involved,  and  could 
serve  as  a  basis  for  more  general  cases  in  which  the  plant  order  may  be  high. 

First,  it  must  be  noted  that  the  quadratic  drag  function  is  in  fact  not  anal3rtic,  so  that 
the  convergence  theorem  above  cannot  apply  formally  to  our  underwater  cable  system.  Al¬ 
though  this  might  lead  to  some  argiunent  over  the  applicability  of  this  approach  to  the 
cable  problem,  one  has  to  remember  that  1)  quadratic  drag  itself  is  invalid  near  the  ori¬ 
gin,  2)  we  will  not  be  carrying  the  infinity  of  terms  required  to  reconstruct  most  analytic 
functions  /(z),  and  3)  we  will  not  be  carrying  out  full-state  feedback.  The  point  is  that 
the  nonlinearity  can  be  approximated  qtiite  well  with  third-order  pol3momials  in  a  given 
regime,  obviously  better  than  a  linear  drag  model.  Our  goeJ  in  Chapter  5  will  be  to  show 
restilts  that  make  it  worthwhile.  For  example,  we  expect  that  situations  requiring  a  very 
large  dynamic  range  will  be  served  well  by  nonlinear  feedback,  as  linearizations  can  be  valid 
only  locally. 

Since  we  are  considering  odd  functions,  the  matrix  Fi  is  zero  and  the  O(z^)  equation  is 

-PiK)P2X^^^  +  P2X^^\A-KPi)z  =  0.  (4.27) 

The  unique  solution  is  Fj  =  0;  this  can  be  verified  numerically  for  any  values  of  x.  Following 
through,  the  O(x^)  equation  is  found  to  be: 


(A^  -  PiK)P3xi^^  -t-  F3arL^’(>l  -  J^Pi)z  +  x^^^^FfPix  +  PjFax^^^  =  0.  (4.28) 

The  third-order  equation  above  can  be  solved  for  all  x  provided  that  we  su^e  able  to 
specify  a  matrix  T[3]  such  that 

z^^\A  -  KPi  )z  =  r[3)zf=’J,  Vr3/?",  (4.29) 

and  another  matrix  ^[3]  such  that 

z^FTp^z  =  Aj3]z'®>,  Vz3F".  (4.30) 
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4.6.S  Definitions  and  Relations  in  Matrix  Calculus 


We  first  need  to  define  a  few  tools  for  matrix  calculus.  A  good  reference  is  Brewer  [20]. 

In  is  an  identity  matrix  of  size  n  x  n.  Let  the  elementary  matrix  be  a  p  X  9  matrix 
with  a  one  in  the  t-th  row,  j-th  column  element,  and  all  other  terms  zero.  Let  A  be  a  p  x  9 
matrix,  and  J?  be  an  r  x  /  matrix.  The  Kronecker  product  [15]  of  A  and  B  is  the  pr  x  ql 
matrix 

oijB  •  •  ■  OlqB 

Repeated  Kronecker  products  are  denoted  by 


A®B  := 


aiiB 

021B 

L 


A*''  :=  A  ®  A  ®  ®  A 

At  times 

The  pq  X  pq  permutation  matrix  Upxq  is 


Example: 


i=i  i=i 


f^Jx3 


1  0  0  0  0  0 
0  0  1  0  0  0 
0  0  0  0  1  0 
0  1  0  0  0  0 
0  0  0  1  0  0 
0  0  0  0  0  1 


The  j-th  column  vector  of  a  matrix  A  is  denoted  as  A,,j.  Then  the  vectored  represen¬ 
tation  of  A  [75],  vcc(A),  is  given  as 
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P.l  B  ^  A  —  Urxp(.A  ®  B)Uqxi 

P.2  {A  ®  B){C  9D)  =  AC9BD 

P.3  vec(AB)  =  (J5^  ®  .fp)t>ec(A) 

P.4  =  + 


Table  4.1:  Properties  of  Kronecker  products  and  matrix  calculus. 


A,, 2 

■^*.9  . 

We  use  the  expression  ‘‘unvec(’)”  to  mean  the  tmstacked  form  of  the  argument. 

The  formal  derivative  of  a  matrix  with  respect  to  a  column  vector  z  of  length  n  is  defined 
as  follows: 


tJ€c(A)  := 


d'A 

dz 


■£Er' 


It  is  important  to  note  that  this  definition  puts  vector  derivatives  of  vectors  at  variance 
with  the  normal  definition,  indeed  with  that  of  Equation  4.28.  The  two  are  related  simply 
by  yj  =  unvec  shall  keep  this  distinction  in  the  sequel.  The  reason  for  not 

using  equivalent  forms  is  that  manipulations  with  the  formal  derivative  are  clearer,  as  will 
be  seen  in  the  sequel. 

Table  4.1  gives  a  number  of  properties  that  will  be  useful  in  our  development.  P.l  is 
due  to  Barnett  [7],  P.2  is  due  to  Bellman  [15],  P.3  is  due  to  Neudecker  [75],  and  P.4  and  P.5 
are  from  Vetter  [100]. 
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4  6.4  Specification  of 

The  O(z^)  vector  in  Equation  4.22  is  taken  to  satisfy  the  implicit  equation 


r®*"  -  1. 

kl 


unvec 


(4.31) 


This  may  seem  a  rather  odd  definition,  and  in  fact  it  is  not  unique.  One  solution  that  works 
is  zIM  =  z®^;  in  this  case,  G  =  (7^,  although  it  is  not  the  identity. 

Example:  Let  n  =  2  and  k  =  2.  Then 


•  • 

2l-«l 

10  0  0 

zW  = 

ZiZj 

ZjZi 

,  and  G  = 

0  M  0 

olio 

Z2Z2 

0  0  0  1 

The  redundancy  in  can  be  eliminated  by  merging  repeated  combinations  in  the 
subscripts  of  the  elements.  We  call  a  thus  reduced  a  minimal  representation  of  z®**. 

n  +  Jb  —  1 


z®*  has  length  n**,  while  zW  has  length  m  := 


[63]. 


We  now  show  that  if  zl^l  is  a  minimal  form  of  z®^,  then  there  is  a  unique  scaling  of  the 
elements  in  z(^I  such  that 


z®*'  =  G^zW  =>  zW  =  (?z®^  (4.32) 

This  fact  is  important  for  our  later  specification  of  T[3].  To  show  that  this  is  true,  we  compute 
GG^zW  explicitly  and  then  substitute.  Let  the  operator  ^/^•...(•)  denote  the  repeated  formal 
derivative  of  the  argument  with  respect  to  Zj,  then  Zj,  and  so  on.  There  are  a  total  of  k 
subscripts.  We  have 


•••  ■ 

■  ,i‘l  ■ 

1 

’  Er=.«n...( ?')*'*' 

' 

1 

^  *! 

• 

J;'). 

(4.33) 
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where  the  ordering  of  derivative  subscripts  is  consistent  with  the  definition  of  the  formal 
derivative.  In  particular,  the  last  subscript  is  changing  slowly  down  the  column  in  the  right- 
hand  side  above.  If  contains  no  repeated  combinations,  then  each  summation  above 
reduces  to  just  one  term.  For  simplification,  let  us  define  an  operator  g[’)  which  takes  an 
integer  into  a  k-tuple  of  integers  defining  the  derivative.  For  example,  $(1)  =  {1, 1,  *  *  *}  and 
g(n^)  =  {n,n,  •  •  •}.  In  addition,  let  the  sorting  operator  s(-)  rearrange  k-tuples  of  integers 
into  ascending  order.  For  instance,  s(l,  1,2,1)  =  {1, 1,1,2}.  Combining  these  observations 
and  definitions  into  Equation  4.33,  we  find  that 


I  rW 


(4.34) 


We  need  two  more  indicial  operators:  let  /(•)  take  a  k-tuple  of  integers  into  an  n-tuple  of 
integers  which  give  the  multiplicity  of  terms  in  the  argument.  For  example,  if  n  =  5  and  k  = 
8,  then  /(1, 1, 2, 2, 1, 3, 4, 2)  =  {3, 3, 1, 1, 0}.  The  number  of  visually  unique  permutations  of 
a  Jb-tuple,  denoted  p(‘),  is  related  to  /  as  follows: 


pi-) 


k\ 

n7=a />(•)!■ 


(4.35) 


Now  we  consider  the  i-th  row  of  (J-  is  nonzero  if  auid  only  if  s{g{j))  =  j(i(2j*^)) 

where  t(-)  is  the  index  Jb-tuple.^  The  number  of  such  occurrences  in  this  row  of  G  is  precisely 
p(i(z|*^)),  and  it  is  clecir  that  these  multiply  with  like-indexed  terms  in  the  vector 
Therefore  the  value  of  the  i-th  row  element  of  GG^zW  is 


(4.36) 


The  derivatives  in  this  equation  are  scalars,  and  if  z|^^  has  associated  with  it  a  coefficient 
Oi,  then 


I, = «.  n  /,(•(?'))! 

^  ^  j=l 


*e.g  ,  =  {1,3,2} 


(4.37) 
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Some  cancellation  then  occtus  in  Equation  4.36,  and  we  find  that 


(4.38) 


Thus,  if  we  set  Oj  =  \/p(i(z|*^)),  then  we  have  the  result  that 


zW  = 


(4.39) 


Simple  substitution  of  Equation  4.31  into  this  form  gives  the  desired  conclusion  that 


zW  =  Gz^^. 


(4.40) 


To  round  things  out,  we  also  have 


z®*  =  G'^Gz^^.  (4.41) 

When  zt*]  is  minimal  and  has  the  proper  scaling.  Equations  4.31,  4.39,  4.40,  and  4.41 
provide  useful  relations  between  z**  and  zW.  We  see  that  Cr^G  may  not  be  the  identity 
because  z®*’  has  repeated  terms.  On  the  other  hand,  GG^  must  be  the  identity  because 
there  are  no  repeated  terms  in  A  special  case  of  our  permissible  z^*’]  is  the  lexicographic 
listing  of  Brockett  [22],  also  used  by  Loparo  and  Blankenship  [63]  and  Yoshida  and  Loparo 
[105].  As  an  example,  for  n  =  3  and  k  =  2,  the  lexicographic  vector  is 


Z^*^  =  [^1  V2ziZ2  V2ziZ3  z|  y/2Z2Z3  zf]^. 


From  this  point  on,  our  developments  incorporate  this  specific  arrangement. 

The  lexicographic  vector  has  a  number  of  other  interesting  properties,  including  the 
following  useful  one: 

Suppose  y  =  Ax,  where  A  is  n  x  n.  Then  the  matrix  equation  y{*l  =  Af*^*f*J  is  satisfied  by 

aW  =  GA®*G^.  (4.42) 


Showing  this  requires  only  a  few  operations:  we  have 
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by  hypothesis, 
=  by  Equation  4.31, 

=  by  P.2, 

=  by  y  =  Az^ 

=  y^*^  by  Equation  4.40. 


4.6.5  Solving  for  T(*] 

Now  we  finally  address  the  problem  of  how  to  pick  a  matrix  Tj*]  such  that 

zWTz  =  TjfcjrW,  (4.43) 

where  T  =■  A  —  KP\  in  the  case  of  Equation  4.29.  An  important  point  is  that 

if  ^  =  Tr,  ,  then  =  rW^z,  (4.44) 


(J„  +  dtT)z{i) 

(/„  +  dtr)WzW(t)=4.  (4.45) 

[(/„  +  dtTp  -  7„]  zW(t)  (4.46) 

Hm  1  [(/„  +  dtTf\  -  /„]  zW  (4.47) 

Is,  i[(^" + 

This  definition  holds  for  arbitrary  n  and  k.  Thus,  we  have  only  to  evaluate  Equation  4.48 
with  T  =  A  -  KP\  and  k  =  3  in  order  to  get  the  matrix  T[i]  of  Equation  4.29.  With  this 
notation,  Equation  4.28  becomes: 

(A^  -  PrK)P3X^^^  +  +  PsTiapW  +  =  q.  (4.49) 


and  following  Brockett  [21],  one  has 

z(t  -I-  dt)  ~ 
zW(t  +  dt)  ~ 
zt*J(t  +  dt)  ~  zW(t)  ~ 

■  ^ 
dt  ~ 

T[k]  := 
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4.6.6  Specification  of  ^(3] 

The  remaining  task  now  is  to  find  a  matrix  such  that  Equation  4.30  is  satisfied.  The 
specification  of  is  somewhat  different  than  that  of  T^s],  since  the  derivative  is 
transposed,  preventing  us  from  using  the  derivative  arguments  above.  A  solution  can  be 
obtained  nonetheless,  and  in  this  section  we  consider  only  the  ib  =  3  case.  As  before,  the 
nature  of  the  reqiiired  calculations  for  more  general  situations  will  be  apparent. 

First,  we  take  three  formal  derivatives  of  Equation  4.30  and  unstack  to  obtain 

-Pi  2]^  =  6A|3]G,  (4.50) 

Now  the  simplifying  substitutions  that  a  :=■  and  6  :=  F^Pi  are  invoked,  and  rising 
Property  P.4,  we  can  expand  out  the  left-hand  side  of  Equation  4.50  as  follows: 


d'^jabz) 

dz^ 


I  (/n  8  (A  8  6)|£  +  (/„  ®  7,  8  |:)  (/.  8  /„  8  6)0  + 

8  6)^  +  (/.  8  (/„  8  /,  8  6)0  +  (4.51) 

d'  d^Z 

^(In  a)(/„  ®  /„  ®  b)j^  -I- 

d'^z 

{In  ®  in  ®  ®  a)(/n  ®  /n  ®  /n  ®  6)-^. 


Since  z  has  only  the  first  derivative  with  respect  to  itself  nonzero,  and  since  a  is  0(z*), 
all  but  the  second,  third,  and  fifth  terms  on  the  right  hand  side  above  are  zero.  Using  the 
properties  P.l  and  P.5  now  to  simplify  Equation  4.50  farther,  we  see  more  nice  symmetry  : 


eXi3]G  = 


(4.52) 


unvec  -J-  7„  ®  £7„x„»  +  /«  ®  ®  f^nxn)  ®  /„  j  t7mxn(/n  ®  b)vec{I„)  . 

There  is  presumably  a  close  relationship  between  ^  ®  /„  and  G;  it  is  a  fact  that  with 
the  current  notation. 
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(4-53) 

This  may  be  verified  for  arbitrary  k  by  simply  writing  out  the  matrix  expressions.  Vf airing 
this  substitution  into  Equation  4.52  above,  and  invoking  property  P.3,  we  have  (in  stacked 
format)  that 


(f^nxn*  +-^n®£^nxn*  + /n  ®  in  ®  t^nXn)(G'^  ®  in)C^mxn(in  ®  6)wec(/„)  =  (G*' 0  in)t»ec(X[3]). 

(4.54) 

The  only  remaining  task  is  to  get  the  two  terms  {G^  0  i„)  to  drop  out.  Since  the  first 
parenthetical  expression  on  the  left  hand  side  of  Equation  4.54  cannot  be  the  identity,  we 
must  show  instead  that  (G^  0  /„)  has  a  left  inverse.  Using  property  P.2  and  the  fact  that 
GG^  =  4n,  it  is  easy  to  see  that  (G  0  /„)  does  the  trick: 

(G  0  i„)(G^  0  in)  =  (GG^)  0  in  =  in.n.  (4.55) 

The  fbal  result  is  then  that 


Vec(A'(3))  =  (G  0  in)(t^nxn>  +  in  ®  f^nXn>  +  4.  ®  in  ®  C^nxn)(<?^  ®  in)Umxn(in  ®  6)vec(i„), 

(4.56) 

which  can  be  solved  using  straightforward  matrix  algebraic  operations. 

A  direct  investigation  of  Equation  4.30  indicates  that  for  k  =  2, 3  at  least,  we  could  also 
use 

n*— n+l  _ 

=  ‘  E  ,  (4.57) 

j=l,n+l,2n+l,— 

where  c^‘*'"“^(’)  means  to  take  only  the  y-th  through  y  +  n  -  1-th  column  of  the  argument. 
Arriving  at  this  conclusion,  however,  requires  a  very  large  amount  of  inspection.  The 
connection  between  Equation  4.57  and  Equation  4.56  can  be  verified. 

Now  the  solution  for  P3  is  at  hand:  Equation  4.28  reduces  to 
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{A^  -  PiK)P3  +  PzT[i)  +  +  PiFz)  =  0.  (4.58) 

This  is  a  standard  equation  in  Sylvester  form. 

4.6.7  Discussion 

This  section  has  presented  a  technique  for  computing  for  arbitrary  plant  order  and 
arbitrary  k.  In  addition,  we  laid  out  the  crucial  machinery  for  finding  under  the  same 
conditions;  even  though  the  solution  was  carried  out  for  h  =  3,  Equation  4.56  is  general 
in  that  when  k  changes,  only  the  summation  of  U's  changes.  An  important  restriction 
of  these  developments  is  that  they  do  not  provide  aU  the  necessary  tools  for  solving  the 
0(z*)  versions  of  Equation  4.25.  An  inspection  of  that  equation  shows  that  some  other 
complicated  terms  arise  even  at  h  =  4;  for  example,  one  of  them  is  Despite  the 

inadequacy  of  the  current  methods  to  treat  such  situations,  it  is  expected  nonetheless  that 
these  problems  will  require  only  extensions  of  our  techniques.  This  woiild  be  an  excellent 
(although  perhaps  rather  tedious)  topic  of  further  research. 

This  method  of  solution  has  the  generality  to  treat  plant  models  of  arbitrary  order.  The 
user  shotild  keep  in  mind  that  although  the  matrix  P3  is  of  size  n  x  m,  some  of  the  inter¬ 
mediate  matrices  in  the  calculations  can  be  extremely  large.  For  example,  Equation  4.56 
requires  the  use  of  square  matrices  containing  n*  elements;  this  would  be  an  rmreasonable 
requirement  for  any  n  above  five  or  so.  Luckily,  the  heuristic  form  of  Equation  4.57  has  no 
such  demands.  The  only  required  matrices  are  G,  which  is  size  m  x  n^,  and  6,  of  size  m  x  n. 
The  largest  remaining  matrix  in  the  algorithm  is  which  has  elements.  Depending 
on  the  hardwue  and  memory  available,  the  plant  order  can  thus  be  relatively  high,  as  will 
be  demonstrated. 

The  matrix  operations  are  relatively  straightforw^lrd;  the  only  nontrivial  operation  is 
solving  the  Sylvester  equation  (4.58).  However,  without  going  into  the  details,  it  has  been 
found  for  the  current  work  that  the  standard  algorithms  are  completely  adequate  (e.g.,  [8]). 

Another  important  design  point  is  that  polynomial  approximations  can  be  valid  only 
in  a  neighborhood  of  the  origin.  Outside  of  this  regime,  the  approximation  and  the  true 
nonlinearity  must  diverge,  which  may  give  rise  to  instability  or  poor  performance  in  the 
closed  loop.  More  will  be  said  about  this  issue  in  Chapter  5. 
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4.7  Design  Example 


The  design  technique  outlined  above  is  ideally  suited  to  plants  with  “hardening”  types  of 
nonlinearity.  A  saturation  nonlinearity,  for  example,  is  difficult  to  express  as  a  low-order 
polynomial  series.  Our  method,  however,  can  accommodate  unstable  plants,  and  here  are 
to  be  found  the  most  dramatic  differences  between  the  nonlinear  optimal  control  solution 
(referred  to  in  the  sequel  as  “NLQR”)  and  the  LQR.  As  a  design  example,  we  consider  a 
plant  based  on  the  first-order  models  of  Chapter  2,  but  incorporate  a  positive  drag  law,  such 
that  near  the  origin  the  drag  forces  are  stabilizing,  but  away  from  it,  they  are  destabilizing. 
The  plant  is  described  by: 
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The  polynomial  approximation  involves  a  linear  plus  a  cubic  part,  and  the  coefficients  for 
these  terms  can  be  found  by  least-squares  minimization.  The  state  penalty  matrix  Q  is 
taken  to  be  C^C  with  C  :=  [0  0  0  0  0  1],  R  =  1,  and  the  initial  displacement  error  is 
one  unit. 

The  LQR  and  NLQR  algorithms  alike  are  sensitive  to  the  linearization  (LQR)  or  domain 
of  polynomial  approximation  (NLQR).  For  this  reason,  we  first  consider  two  designs  with 
no  model  error,  based  on  two  speed  regimes.  In  Figure  4-3  is  shown  the  simulated  response 
of  the  closed-loop  system  for  speed  envelopes  of  0.7  (top)  and  3.5  (bottom).  The  value 
of  0.7  is  the  maximum  speed  reached  in  the  NLQR  response;  it  is  evident  that  when  the 
envelope  is  minimized  in  this  way,  the  LQR  is  the  superior  design.  The  cost  integrals  over 
twenty  seconds  are  J  =  3.21  and  J  =  5.60  for  the  LQR  and  NLQR  responses,  respectively. 
This  may  be  regarded  as  a  very  good  linearization  from  the  point  of  view  of  the  LQR,  and 
a  tenuous  one  for  the  NLQR,  as  any  plant  speeds  outside  the  envelope  access  an  inaccurate 
cubic  drag  in  the  control  law. 
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One  may  extend  the  speed  envelope  to  encompass  all  the  motions  that  might  be  encoun¬ 
tered,  giving  the  NLQR  a  greater  operating  range.  The  effects  of  this  design  modification 
applied  to  the  same  plant  are  illustrated  in  the  lower  half  of  Figure  4-3.  We  observe  that 
the  NLQR  response  is  largely  the  same  as  in  the  previous  case  (7  =  4.59),  but  that  the 
LQR  suffers  drastically,  as  linear  damping  cannot  apply  to  the  entire  speed  envelope.  It  is 
interesting,  however,  that  rather  than  exhibiting  instability,  the  LQR  limits  the  response 
speeds  to  well  within  the  regime  of  negative  drag. 

The  issue  of  robustness  provides  another  scenario  for  comparison  of  the  LQR  and  NLQR. 
For  example,  the  runs  of  Figure  4-3  have  been  repeated  in  Figure  4-4,  with  the  following 
modeling  errors  incorporated  in  the  design  plant:  1)  underestimation  of  mass  by  one-half,  2) 
overestimation  of  stiffness  by  a  factor  of  two,  and  3)  imderestimation  of  the  “drag  coefficient” 
by  one-half,  to  be  applied  as  i(|  r,  |  —1).  In  the  case  of  the  0.7  speed  envelope,  the  LQR 
loop  is  unstable,  while  the  NLQR  leads  to  approximately  the  same  result  before  (7  = 
2.77).  The  NLQR  performs  well  in  the  3.5  speed  envelope  case  as  well  (7  =  1.98),  and  the 
LQR  repeats  the  slow  but  stable  response  of  Figure  4-3. 

The  effects  of  different  envelope  sizes,  modeling  uncertainties,  and  other  topics  are  dis¬ 
cussed  further  in  Chapter  5,  using  experimental  data.  For  the  example  of  this  section,  it  is 
apparent  that  the  NLQR  has  a  significcintly  larger  range  of  operation  than  the  LQR,  and 
has  improved  robustness  properties,  as  expected. 


4.8  Proposed  Synthesis  Techniques 

We  have  developed  the  necessary  machinery  for  model-based,  closed-loop  control  of  finite¬ 
dimensional  models  of  a  towed  cable  system.  The  approaches  are  bxiilt  upon  the  LOR 
technique  for  nonlinear  control,  and  thus  benefit  from  nonlinear  estimators  and  control 
laws.  The  complete  controller  can  be  synthesized  into  numerous  combinations;  six  of  them 
are  proposed  for  the  towing  problem: 

1.  Linear  loopshaping. 

2.  LTR:  LQR  with  KF 

3.  LOR:  LQR  with  CGEKF. 

4.  LOR:  LQR  with  EKF. 
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displacement  displacement 


Figtire  4-4:  Simulations  of  LQR  and  NLQR  responses  for  speed  envelopes  of  0.7  (top)  and 
3.5  (bottom),  with  modeling  imcertainty. 
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5.  LOR:  Optimal  nonlinear  control  with  CGEKF. 

6.  LOR:  Optimal  nonlinear  control  with  EKF. 

Linear  loopshaping  (e.g.  [38])  is  to  he  considered  because  it  is  by  far  the  easiest  com¬ 
pensator  to  design  for  a  SISO  linear  plant  model.  The  LTR  method  can  he  expected  to 
provide  control  nearly  identical  to  that  of  loopshaping,  but  is  computed  nonetheless  since 
its  parameters  are  used  directly  in  the  nonlinear  components  which  follow.  The  proposed 
nonlinear  designs  cover  the  range  of  possible  combinations  valid  in  LOR,  and  are  listed  in 
order  of  increasing  computational  requirements.  As  previously  noted,  the  nonlinear  optimal 
feedback  law  is  denoted  “NLQR”  in  the  sequel. 
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Chapter  5 


Closed-Loop  Experiments 


5.1  Introduction 

In  this  chapter,  we  apply  the  theory  of  Chapters  3  and  4  in  order  to  physicaUy  control  scale 
models  of  two  relevant  physical  systems.  Since  the  plant  is  known  to  be  naturally  stable, 
the  main  emphasis  is  not  on  stability,  but  on  the  myriad  of  other  properties  to  which  the 
designs  aspire.  These  include  transient  performance  and  long-term  behavior,  robustness, 
disturbance  rejection,  and  many  other  points. 

The  physical  systems,  scaling  procedure,  and  experimental  apparatus  are  described  first. 
Then  some  practical  implementation  issues  are  discussed,  such  as  model  identification  and 
controller  design  and  implementation.  Experimental  resiilts  in  regulation  and  closed-loop 
tracking  are  provided  for  the  scaled  systems,  and  the  general  trends  are  discussed. 

5.2  Experimental  Setup 

5.2.1  Specification  of  Two  Nominal  Full-Scale  Systems 

The  model  tests  were  scaled  to  study  the  behavior  of  two  full-scale  systems  of  interest. 
The  first  is  a  version  of  the  ARGO/JASON  system,  which  has  a  very  slender  cable  and  a 
significant  clump  weight.  The  second  is  a  drilling  deployment,  which  entails  a  very  large, 
stiff  pipe,  and  a  comparatively  small  towfish  at  the  lower  endpoint.  The  physical  parameters 
for  these  two  nominal  systems  are  given  in  Table  5.1. 
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ARGO/JASON 

Drill  String 

Vessel  maximum  speed 

1  m/s  towing 

1  m/s  towing 

Vessel  bandwidth 

.02  rad/ s 

.02  rad/» 

Cable  diameter 

1.73  cm 

41  cm 

Cable  eff.  mass 

1.18  kg/m 

378  kg/m 

Cable  weight  in  water 

6.9  N/m 

1118  N/m 

Cable  drag  coefficient 

1.6 

1.6 

Cable  El 

878  N  •  m*i 

76e6  N  ■  m* 

Depth 

1000  m 

1000  m 

Towfish  effective  mass 

900  kg 

1000  kg 

Towfish  weight 

6900  N 

mo  N 

Towfish  area 

1  meter^ 

2  meter^ 

Towfish  drag  coefficient 

2 

2 

Table  5.1:  Physical  parameters  for  two  nominal  full-scale  systems. 

5.2.2  Laboratory  Scale-Models 

The  nondimensionalization  of  the  first-order  lateral  cable  equations  in  Section  2.5.3  leads 

directly  to  a  set  of  nondimensional  par  [  PrinterError:  out  of  paper  )  ameters  that  must  be  maintair 

also  [79]).  They  are  repeated  here: 


*1 

^2 

*4 

6 

e 


(m+ma)g' 

W 

ma)g' 

W 


{M  +  Ma)g' 
pCj(P 
2(m  +  TTla)  ’ 
pCdtAjd 

2{M  +  M^y 


and 


It  is  noteworthy  that  the  variable  c,  which  was  not  included  in  Chapter  2,  is  the  ratio  of 
nonlinear  vehicle  drag  forces  to  vehicle  inertial  forces.  In  addition  to  these  five  parameters, 
the  Reynolds  number  should  be  maintained  through  scaling,  and  the  amplitude-to-diameter 
ratio  should  be  kept  large.  These  steps  ensure  that  we  are  not  working  with  two  drastically 
different  flow  regimes.  In  addition,  the  bending  stiffness  parameter  of  Section  2.2.3  should 
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be  small  for  both  the  full-size  and  the  scaled  systems. 

A  test  rig  was  built  at  the  Coastal  Research  Laboratory  of  the  Woods  Hole  Oceano¬ 
graphic  Institution  to  carry  out  the  scale-model  tests.  The  physical  arrangement  is  shown 
in  Figure  5-1.  An  I-beam  across  the  surface  of  the  3.3-meter-deep  test  tank  was  fitted  with 
a  rolling  carriage,  which  served  as  the  upper  endpoint  for  the  cable  in  the  water,  and  was 
positioned  through  a  chain  drive  and  motor,  controlled  from  a  personal  computer.  The 
carriage  position  was  tracked  using  the  resolver  in  the  motor,  and  the  lower  endpoint  of  the 
cable  was  fitted  with  a  small  acoustic  transducer.  This  transducer  communicated  with  a 
short-baseline  acoustic  net  in  the  300  kHz  range,  mounted  in  the  test  tank. 

The  first  experimental  setup  was  designed  to  approximate  the  ARGO/JASON  fiiU-scale 
system.  For  the  cable,  the  coaxial  link  from  the  transducer  to  the  surface  was  adequate,  and 
the  transducer  alone  gave  a  good  approximation  of  the  ARGO/JASON  towfish.  Table  5.2 
shows  how  well  the  scaling  held.  As  required,  the  bending  stifihess  was  kept  small,  the 
Reynolds  mamber  was  in  the  same  flow  regime  as  the  full-scale  system,  and  the  amplitude- 
to-diameter  ratio  was  high. 

The  second  experimental  setup  was  a  scaled  version  of  the  full-scale  drill-string  system. 
In  this  case,  the  much  larger  diameter  and  weight  required  the  use  of  half-inch  Tygon 
tubing,  with  two  3.2-millimeter  wire  ropes  and  the  coaxial  cable  inside  for  mass  and  weight. 
In  order  to  eliminate  the  significant  in- water  weight  of  the  transceiver  at  the  bottom  (0.57 
Newtons),  foam  flotation  was  added  at  the  bottom  of  the  cable.  This  foam  introduced 
significant  frontal  area  and  added  mass  to  the  transducer,  but  made  the  overall  scaling  in 
the  cable  itself  much  better.  In  other  words,  we  could  achieve  good  scaling  in  the  cable  or 
at  the  endpoint,  but  not  both;  the  former  was  chosen  for  this  work.  In  Table  5.2,  we  see 
that  the  parameters  k4  and  c  are  in  significant  error. 

Another  important  note  is  that  flow  aroxmd  a  full-scale  drill  string  may  be  supercritical 
at  one  meter  per  second,  while  in  the  test  tank,  we  were  \mable  to  match  this  regime  because 
of  velocity  constraints.  In  supercritical  flow,  such  a  pipe  could  have  a  substantially  reduced 
drag  coefficient,  especially  if  its  surface  is  smooth  [81].  This  scaling  disparity  may  or  may 
not  be  significant  from  a  practical  control  systems  point  of  view,  and  it  was  decided  to  use 
the  scale  model  as  described  above,  given  the  difficulty  of  matching  all  the  nondimensional 
parameters.  In  the  experiments  which  follow,  however,  we  do  consider  the  case  where  the 
modeled  drag  terms  are  overestimated  by  a  factor  of  three. 
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I-beam  with  cable  drive 


3.3-meter  tank 


Figure  5-1:  Schematic  view  of  the  WHOI  test  tank. 


ARGO/JASON 

Scale  A/J 

Drill  String 

Scale  DS 

J.  de  Fuca 

.59 

.59 

.27 

.64 

*2 

.60 

.79 

.0022 

.0017 

.21 

*4 

.78 

.73 

.82 

.014 

.75 

€ 

.21 

.21 

.36 

.37 

.16 

£ 

.020 

.027 

.83 

.22 

.0022 

El 

iv,^L+W)Li 

7.5e-9 

4.0e-5 

2.1e-5 

2.2e-4 

l.le-10 

Re 

1.3e4 

1.8e3 

3.1e5 

9.8e3 

4.9e3 

Aid 

5-8€3 

4.2e2 

2.5e2 

7.9el 

1.6e4 

Table  5.2:  Scaling  for  model  tests. 

With  regard  to  the  lateral  motions  themselves,  since  the  ratio  of  excitation  to  diameter 
was  small  for  both  the  full-scale  and  model  systems,  the  important  nondimensional  variable 
became  the  angle  6.  As  ^  was  at  least  an  order  of  magnitude  smaller  in  the  test  tank  than  in 
the  full-scale  systems,  preserving  6  meant  that  the  angle  had  to  be  that  much  larger  in  the 
test  tank.  Thus,  a  one-meter  excursion  in  the  test  tank  corresponds  to  approximately  thirty- 
three  meters  for  the  full-scale  drill  string,  and  just  seven  for  the  ARGO/JASON  system.  In 
addition,  the  full-scale  frequencies  correspond  with  model  frequencies  approximately  twenty 
times  as  fast.  Thus,  a  vessel  with  a  DP  bandwidth  of  0.05  radians  per  second  implies  an 
upper  frequency  limit  of  one  radian  per  second  in  the  test-tank. 

Vertical  excursions  of  the  towfish  were  not  scaled  accurately  in  our  test  tank,  since  we 
only  scaled  the  first-order  cable  equations,  which  entail  no  vertical  coupling. 

The  full  listing  of  physical  parameters  for  the  scale  experiments  of  this  thesis  is  given  in 
Table  5.3. 

5.2.3  Dividing  the  Systems 

Table  5.2  indicates  a  significant  physical  difference  between  an  ARGO/JASON-typc  deploy¬ 
ment  and  a  drill  string  deployment,  in  particular,  the  fraction  which  is  the  the  ratio  of 
towfish  weight  to  the  weight  of  the  cable,  differs  by  two  orders  of  magnitude.  In  practical 
terms,  this  means  that  the  cable  shape  for  a  towed  ROV  system  is  dominated  by  the  weight 
of  the  towfish,  and  the  cable  is  generally  concave  downwards.  For  a  drill  string,  the  lower 
portions  of  the  pipe  have  decreasing  resistance  to  drag  forces  because  of  the  small  tension; 
the  end  of  the  pipe  is  comparatively  free  to  “whip.”  As  a  result,  the  shape  of  a  towed  drill 
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Scale  DS 


Scale  A/J 

Scale  DS 

Maximum  towing  speed 

1  m/s 

<= 

Maximum  excursion 

3.5  m 

Cart  tracking 

motor  resolver 

(o'  =:  .1  mm) 

Actuation  bandwidth 

3  rad/s 

<= 

Cable  material 

2.4  mm  coax 

12.7  mm  Tygon 
with  steel  wire 

Cable  effective  mass 

.022  kg/m 

.353 

Cable  weight /length  in  water 

.13  N/m 

.92 

Nominal  cable  C4 

1.6 

1.6 

Depth 

2.74  m 

Towfish 

transceiver 
.57  N  weight 

transceiver 
and  foam,  .05  N 

Towfish  effective  mass 

.080  kg 

.15 

Towfish  area 

9  cm* 

25 

Nominal  towfish  Cj 

2 

2 

Towfish  tracking 

acoustic 

(0  ~  1  mm) 

Sample  rate  (no  computations) 

0.22  seconds 

•<= 

Sample  rate  (computations) 

0.275  seconds 

<= 

Pure  delay 

^  0.5  seconds 

-  1.6 

Table  5.3:  Physical  parameters  for  WHOI  test  tank  experiments. 
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string  is  usually  concave  upwards. 

These  two  systems  differ  not  only  in  their  physical  structure,  but  also  in  their  usual 
operation.  A  towfish  which  is  part  of  an  ROV  system  often  has  a  tether  of  length  fifty  or 
more  meters  connecting  it  to  the  ROV,  so  high  accuracy  positioning  is  not  crucial.  As  a 
sensor  platform,  however,  the  towfish  may  be  required  to  follow  predetermined  tracklines. 
The  drill  string  presents  a  different  operational  need:  there  is  usually  no  requirement  for 
trajectory-following,  but  regulation  of  the  endpoint  position  is  crucial. 

The  fact  that  these  two  full-scale  plants  are  so  different  in  terms  of  their  physical  response 
suggested  experimental  tests  with  two  scaled  systems.  For  example,  it  was  unclear  what 
the  effects  of  endpoint  whipping  would  be  on  closed-loop  stability.  At  the  same  time,  the 
operational  differences  pointed  to  a  natural  division  between  the  types  of  experiments  to  be 
conducted.  In  the  data  of  this  chapter,  we  therefore  present  regulation  runs  for  the  scaled 
drill  string  system,  and  trajectory-following  runs  for  the  scaled  ROV  system.  This  division 
limits  the  results  for  each  system,  although  it  has  been  verified  that  the  drill  string  motions 
can  be  successfully  preshaped,  and  that  regulation  with  the  ROV  system  can  be  very  good 
also. 


5.3  Practical  Issues 

5.3.1  Drag  Coupling  and  Polynomial  Approximation 

The  restriction  of  cable  motions  to  a  single  plane  is  unrealistic  because  physical  distur¬ 
bances  act  in  all  directions,  and  because  we  may  want  to  follow  reference  trajectories  with 
turns  or  curves.  The  first-order  three-dimensional  equations  are  coupled  only  throiigh  drag, 
and  the  drag  force  coupling  of  Equation  2.38  indicates  that  motions  in  one  plane  tend  to 
increase  the  effective  drag  seen  in  the  other  plane.  For  this  reason,  one  would  not  expect 
dramatic  changes  in  performance  if  the  in-plane  and  out-of-plane  control  problems  were 
addressed  separately  using  the  techniques  discussed  so  far.  Indeed,  if  linear  techniques  are 
used,  then  one  sees  this  procedure  immediately,  since  a  linear  drag  law  has  no  coupling. 
Technically,  this  makes  the  problem  multi-input /multi-output,  although  it  is  clear  that  the 
two  subsystems  are  completely  separate,  and  SISO  loopshaping  can  be  used. 

Incorporation  of  the  drag  coupling  is  not  difficult  in  the  context  of  our  nonlinear  ap¬ 
proaches.  In  Safonov’s  theorems  for  nominal  nondivergence  of  the  CGEKF  (Equation  4.9) 
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and  stability  of  the  LQR  (4.21),  the  drag  coupling  leads  to  slightly  different  tolerances, 
but  the  bounds  are  roughly  the  same  as  in  the  single-plane  case.  The  EKF  mountains  its 
robustness  properties  as  well,  so  long  as  the  estimator  model  contains  the  coupled  forces. 

The  drag  coupling  may  be  easily  included  in  the  nonlinear  optimal  approximation  of 
Section  4.6.  The  key  is  to  modify  the  formulation  for  so  that  the  appropriate  coupling 
is  achieved;  the  rest  follows  easily  within  the  framework  of  the  chapter.  Toward  this  end, 
perhaps  the  easiest  way  to  generate  the  coefficients  is  to  carry  out  a  least  squares  fit  for  grid 
points  on  a  subset  of  the  z  —  y-plane.  For  example,  for  points  over  a  disc  of  unit  radius, 
one  obtains: 


z^z^  +  y*  ~  0.34132  -i-  0.6896**  4-  0.6865zy*, 

leading  to  a  maumum  difference  between  the  functions  of  0.033,  with  an  RMS  error  of 
0.012.  In  the  best-fit  linear  approximation,  an  2-coefficient  of  0.7963  gives  a  maximTun 
error  of  0.19  and  an  RMS  error  of  0.081.  For  the  in-plane  case,  we  have 

2  1  2  |~  0.30402  -I-  0.74932*, 

with  a  maximum  error  of  0.029  (0.018  RMS).  The  coefficient  in  the  linear  case  is  0.7308; 
the  maximum  error  is  0.21,  with  an  RMS  value  of  0.11.  In  general,  it  is  apparent  that  the 
higher-order  approximation  improves  the  overall  accuracy  by  a  factor  of  six,  roughly. 

5.3.2  System  Identification 

The  two  scaled  systems  were  modeled  by  minimizing  the  open-loop  simulation  error  [55], 
as  discussed  in  Section  2.6.  The  input  data  set  used  to  create  the  I/O  data  was  generated 
through  the  following  technique: 

1.  A  random  series  with  a  Gaussian  probability  density,  zero  mean,  and  unity  variance 
was  padded  with  zeros  at  the  start  and  the  end. 

2.  The  series  was  filtered  using  a  second-order  digital  Butterworth  filter,  with  a  cutoff 
frequency  at  three  radians  per  second.  This  frequency  scales  to  0.15  radians  per 
second,  and  was  well  within  the  bandwidth  of  the  apparatus. 
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3.  An  absolute  multiplier  was  applied  to  the  series  so  that  the  full  dimensions  of  the  tank 
were  utilized. 

4.  The  TtiAxinniTn  velocity  was  checked  so  as  to  not  exceed  the  capabilities  of  the  appa¬ 
ratus.  This  speed,  about  one  meter  per  second,  corresponds  to  approximately  1.65 
meters  per  second  for  the  full-scale  drill  string  system,  and  0.35  meters  per  second  for 
the  full-scale  ROV  system. 

The  weight  terms  were  calculated  beforehand,  leaving  the  cable  and  endpoint  drag,  and 
the  cable  and  endpoint  effective  mass  as  unknowns  to  be  identified;  then  system  identifi¬ 
cation  was  carried  out  as  described  in  Chapter  2.  The  performance  of  two-  and  four-node 
models  for  the  scaled  drill  string  system  are  shown  in  Figure  5-2;  the  four-node  model  is 
clearly  superior. 

5.3.3  Quantification  of  Modeling  Errors 

The  models  generated  by  the  above  method  come  with  no  quantitative  measures  of  goodness 
that  are  useful  for  control.  Since  the  robustness  properties  for  optimal  controller  loops  are 
known  (see  Section  4.5.2),  we  now  consider  the  problem  of  quantifying  the  imcertainty 
operator  A.  In  particular,  the  fundamental  measme  of  the  uncertainty  is  the  L2-%^  of  A. 

The  gain  of  a  nonlinear  operator  usually  must  be  found  from  actual  response  (simulation 
or  physically  obtained)  data,  although  clearly  no  amount  of  experimenting  can  guarantee 
that  the  maximum  value  has  been  found.  This  is  in  contrast  to  a  linear  operator  G(b>),  for 
which  the  fr2*i^orm  of  the  output  is  related  to  the  X2'norm  of  the  input  by 
II  G(o>)  llooll  u(f)  II2  [38].  It  has  been  typical  in  the  literature  to  estimate  the  gain  of 
nonlinear  operators  by  forcing  them  with  swept-sine  inputs  of  various  amplitudes.  While 
this  approach  obviously  does  not  yield  much  information  about  the  gain  during  transients, 
like  the  describing  function  technique  with  harmonic  excitation  it  aspires  to  cover  all  of  the 
regimes  through  which  an  unstable,  oscillatory  trajectory  co\ild  pass.  For  this  reason,  we 
consider  only  swept-sine  wave  signals  as  inputs  in  this  section. 

Another  difficulty  is  that  the  operator  we  are  interested  in  is  located  at  the  plant  input, 
consistent  with  our  recovery  of  the  controller  loop.  As  a  result,  we  cannot  compare  the 
output  of  our  model  with  the  output  of  the  actual  plant  for  a  given  input  and  get  useful 
numbers.  We  require  a  completely  new  computational  approach.  Inspection  shows  that  the 
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Figure  5-2:  Identification  of  two-node  (top)  and  four-node  (bottom)  models  for  the  scale 
drill  string  system. 
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uncertainty  operator  A  (see  Figure  4-2)  can  be  expressed  in  terms  of  the  actual  plant  Ga, 
and  the  nominal  plant  for  pre-multiplicative  and  pre-division  uncertainty  structures,  as 
follows: 


G~^Ga  -  I  •  pre-multiplicative 

(5.1) 

G~^Gn  —  I  •  pre-division. 

(5.2) 

It  is  apparent  that  G~^  can  be  calculated  using  our  dynamic  inversion  scheme  of  Chapter  3, 
and  Ga  can  be  found  by  high-order  simulation  or  actual  experiments.  This  immediately 
gives  us  a  way  to  quantify  ||  A  ||2,r  explicitly,  for  a  pre-multiplicative  structure.  A  procedure 
for  calculating  it  is  outlined  as  follows. 

1.  It  was  shown  in  Section  2.5.4  that  the  plant  gain  increases  monotonically  as  the 
excitation  amplitude  decreases,  and  might  become  quite  large  for  small  excitations.  || 
A  ||2,r  also  might  be  quite  large  in  this  regime,  so  a  lower  limit  of  excitation  amplitude 
should  be  picked,  below  which  the  plant  response  is  inconsequential. 

2.  Generate  an  input  signal  u  that  satisfies  the  lower  bound  and  covers  the  relevant 
frequency  range.  A  swept-sine  signal  with  sections  at  different  amplitudes  would  be 
an  effective  choice. 

3.  Drive  the  plant  or  a  high-order  model  (Ga)  with  this  u.  The  output  is  y  :=  GaU. 

4.  Solve  the  inverse  dynamics  for  y,  yielding  u  :=  G~^GaU.  Thus  u  —  u  =  Au. 

5.  Define  a  “notching”  operator  nu(u)  whose  support  is  confined  to  a  single  period  of 
oscillation,  the  frequency  being  u>.  Thus,  ut  :=  ^^(u)  is  a  truncated  version  of  u, 
whose  tails  outside  a  single  period  are  set  to  zero.  Similarly,  let  Ut  be  a  version  of  u 
with  the  same  tails  set  to  zero.  Then  we  have  the  estimate  ut  —  tq  Aut,  and 


II  A  ||2.T=i  max 


II  fo  -  “t  Hz.r 
IKIkr  • 


(5.3) 


Clearly,  the  above  calculation  is  accurate  only  if  the  rate  of  change  of  frequency  in  the 
swept  sinusoid  is  slow  enough  to  accommodate  nearly  steady-state  plant  responses. 
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Figure  5-3:  Verification  of  dynamic  inversion  scheme  for  calculation  of  A. 

The  notching  operator  is  useful  because  the  gains  are  computed  for  complete  periods 
of  oscillation. 

First,  we  give  an  illustration  that  the  inversion  scheme  is  up  to  the  task  by  letting 
Gn  =  Ga  for  the  scaled  drill-string  system.  As  outlined  above,  the  swept-sine  input  u  was 
run  through  Ga,  then  y  through  .  Figure  5-3  shows  u,  v,  and  y;  the  plant  has  clearly 
been  inverted  properly,  as  the  difference  between  u  and  u  is  almost  imperceptible. 

Swept-sine  inputs  and  experimental  responses  for  the  scaled  drill  string  system  are  given 
for  full-scale  amplitudes  of  one,  two,  five,  and  ten  meters,  for  full-scale  frequency  ranges 
between  .01  and  .1  radians  per  second,  in  Figure  5-4.  An  example  inversion  for  the  third 
graph  of  Figtire  5-4  is  shown  in  Figure  5-5,  demonstrating  the  usual  deterioration  of  model 
acctiracy  as  the  frequency  increases.  It  is  intuitive  that  small  variations  in  the  experimental 
towfish  response  can  lead  to  large  variations  in  v;  since  the  forward  system  is  a  low-pass 
filter,  we  expect  the  inverse  system  to  amplify  high-frequency  inputs  to  some  degree. 
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Figure  5-5:  Example  of  u  and  v  for  experimental 
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Figure  5-6:  Gains  of  uncertainty  operator  as  a  function  of  frequency. 


Once  the  inversions  are  completed  for  the  range  of  amplitudes,  the  Lj  gain  of  A  may  be 
computed  according  to  Equation  5.3,  for  the  various  frequency  regimes.  These  frequency- 
dependent  gains  are  plotted  for  the  scale  drill  string  system,  in  Figure  5-6.  Since  by  Sec¬ 
tion  4.5.2,  closed-loop  stability  of  the  nonlinear  optimal  control  loop  is  guaranteed  when 
II  ^  lU^  21  VC  see  that  the  modeling  errors  are  acceptable  below  an  excitation  frequency  of 
approximately  0.5  radians  per  second. 

5.3.4  Controller  Design 

The  crossover  frequency  for  the  nominal  linear  design  should  be  selected  in  a  manner  consis¬ 
tent  with  the  actuation  limits,  and  so  that  the  basic  guidelines  of  loopshaping  and  LTR  are 
followed.  Typically,  these  pertain  to  a  smooth  loop  transfer  function  with  reasonable  phase 
loss  in  the  crossover  region,  and  a  stable  compensator.  For  the  scale  drill  string  system, 
the  maximum  crossover  is  approximately  0.2  radians  per  second,  which  compares  well  with 
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the  modeling  error  analysis  above,  and  the  first  natural  mode  of  the  scaled  plant  located  at 
about  0.8  radians  per  second. 

For  a  given  plant  model,  and  speed  envelope,  an  entire  family  of  compensators  can  be 
constructed  as  follows.  The  loopshaping  and  LTR  controllers  are  first  designed  to  the  same 
specifications,  and  the  nonlinear  designs  are  obtained  through  straightforward  extensions; 

1.  The  gain  matrix  H  from  the  KF  is  used  in  the  CGEKF. 

2.  The  process  noise  gain  matrix  L  for  the  EKF  is  set  to  B,  and  the  sensor  noise  /x  is 
the  same  as  that  used  in  the  LTR  design. 

3.  The  nonlinear  optimal  controller  utilizes  the  same  value  of  p  as  in  the  LTR  design, 
and  has  Q  =  C^C. 

Thus,  the  LTR  forms  the  link  between  the  basic  characteristics  of  the  loopshaping  design, 
and  the  parameters  for  the  nonlinear  components  in  the  more  complex  controUers.  An 
example  LTR  design  loop  transfer  function  is  shown  in  Figure  5-7;  the  linearization  speed 
for  this  case  is  0.2  meters  per  second. 

5.3.5  Formal  Properties  of  the  CGEKF  and  LQR 

In  Chapter  4,  it  was  noted  that  the  CGEKF  and  LQR  can  be  applied  to  some  nonlinear 
systems  with  guaranteed  results  [82].  In  general,  we  must  satisfy  a  bound  on  the  size  of  the 
discrepancy  between  A  for  the  linearized  system,  and  V/  for  the  nonlinear  system.  The 
conditions  are  to  be  checked  through  the  operating  regime.  Figure  5-8  gives  an  indication 
of  how  these  hold  for  the  scaled  drill  string  system  for  the  plant  model  and  controller  design 
of  Figure  5-7,  in  the  presence  of  quadratic  drag,  with  a  linearization  speed  of  0.1  meters  per 
second.  For  each  speed  value  on  the  abscissa,  actual  nodal  speeds  are  taken  as  uniformly 
random  nvimbers  between  fifty  percent  and  one-hundred-fifty  percent  of  the  speed  on  the 
axis.  Overall,  we  see  that  the  conditions  are  generally  met  above  the  linearization  speed, 
but  not  below  it.  Hence,  the  only  way  to  obtain  guaranteed  properties  from  the  CGEKF 
and  LQR  with  respect  to  the  drag  nonlinearity  is  to  set  the  linearization  speed  extremely 
low.  This  requirement  usually  has  the  effect  of  slowing  down  the  transient  response  of  the 
closed-loop  system,  as  will  be  verified. 
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Figure  5-8:  Formal  guarantees  for  LQR  (top)  and  CGEKF  (bottom)  for  a  design  lineariza¬ 
tion  of  0.1  meters  per  second. 


0.1  m/s  ramp  0.2  m/s  ramp  0.4  m/s  ramp 

0.1  m/s  envelope 
0.2  m/s  envelope 
0.4  m/s  envelope 

Figure  5-9  Figure  5-10  Figure  5-11 

Figure  5-12  Figure  5-13  Figure  5-14 

Figure  5-15  Figure  5-16  Figure  5-17 

Table  5.4:  Combinations  of  ramp  speed  and  design  envelope  for  experimental  step  responses. 


5.4  Step  Responses 

In  this  section,  the  ramped  setpoint  change  is  viewed  as  a  generic  operation  for  the  drill 
string  system,  verifying  the  stability  of  the  closed-loop  system,  and  revealing  the  transient 
response  and  long-term  behavior  (i.e.,  limit  cycling).  We  consider  changes  in  the  setpoint  of 
2.0  meters,  corresponding  to  66-meter  moves  for  the  full-scale  drill  string  system.  All  of  the 
designs  have  a  crossover  frequency  of  0.2  radians  per  second,  and  are  based  on  the  two-node 
model.  In  order  to  establish  the  basic  tradeoffs  for  the  different  design  approaches,  nine 
combinations  of  speed  envelope  and  ramp  speed  for  the  setpoint  are  considered;  these  are 
listed  in  Table  5.4. 

The  responses  in  each  figure  are  arranged  in  order  of  increasing  complexity  and  com¬ 
putational  load  for  the  control  and  estimation  algorithms.  Each  experiment  initiates  with 
a  ten-second  dead  time,  to  allow  settling  of  the  physical  system,  and  a  ten-second  period 
in  which  the  controller  is  enabled  but  the  setpoint  is  still  at  the  origin.  The  duration  of 
the  runs  is  chosen  to  be  long  enough  that  the  limit  cycling,  if  any,  is  completely  developed. 
Finally,  as  stated  earlier,  the  controllers  all  have  integral  action,  so  no  steady-state  offsets 
are  expected.  However,  as  the  controllers  are  only  of  Type  1  (a  single  integrator),  tracking 
errors  during  steady  towing  are  not  corrected. 

There  are  a  number  of  general  observations  which  can  be  made  about  the  designs, 
based  on  these  data.  The  most  striking  point  is  that  the  EKF  in  every  case  leads  to 
reduced  limit-cycling  behavior.  In  contrast,  it  is  known  that  the  linear  Kalman  filter  has  no 
nondivergence  properties  when  applied  to  nonlinear  systems,  and  that  the  CGEKF  is  valid 
only  for  operating  speeds  above  the  linearization.  These  facts  are  corroborated  by  the  data, 
and  suggest  that  all  other  results  aside,  the  EKF  is  a  crucial  part  of  the  best  controller 
designs,  with  LOR. 

It  is  apparent  that  since  the  physical  plant  is  very  stable,  one  can  achieve  asymptotically 
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Figure  5-9:  Step  responses  with  0.1  m/s  envelope,  0.1  m/s  ramp  speed. 


Figure  5-10:  Step  responses  with  0.1  m/s  envelope,  0.2  m/s  ramp  speed. 
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Figtire  5-11:  Step  responses  with  0.1  m/s  envelope,  0.4  m/s  ramp  speed. 


Figure  5-12:  Step  responses  with  0.2  m/s  envelope,  0.1  m/s  ramp  speed. 
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guxe  5-13:  Step  responses  with  0.2  m/s  envelope,  0.2  m/s  ramp  speed. 


Figure  5-14:  Step  responses  with  0.2  m/s  envelope,  0.4  m/s  ramp  speed. 
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Figure  5-15:  Step  responses  with  0.4  m/s  envelope,  0.1  m/s  ramp  speed. 


Figure  5-16:  Step  responses  with  0.4  m/s  envelope,  0.2  m/s  ramp  speed. 
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Figiire  5*17:  Step  responses  with  0.4  m/s  envelc^e,  0.4  m/s  ramp  speed. 

stable  results  from  any  linear  design  technique  by  making  the  design  envelope  very  small 
(e.g.,  Figixre  5-9).  We  recall  that  Safonov’s  conditions  validate  both  the  LQR  and  the 
CGEKF  in  this  situation,  and  it  seems  that  the  Kalman  filter  works  well  also.  As  a  general 
approach,  however,  the  specification  of  a  small  envelope  implies  poor  performance  far  away; 
that  is,  at  high  speeds.  This  trend  may  be  seen  in  Figure  5-11,  where  the  loopshaping  and 
LTR  responses  contain  a  substantial  overshoot. 

The  NLQR  presents  a  special  problem  because  when  its  design  envelope  is  exceeded,  the 
cubic  term  in  the  design  plant  model  rapidly  overtakes  the  actual  drag  in  the  plant.  Thus, 
the  plant  has  substantially  less  drag  than  the  control  law  expects.  This  fact  is  in  contrast 
to  the  LQR,  in  which  the  plant  has  more  drag  than  the  control  law  expects  outside  the 
envelope.  Unlike  the  general  deterioration  of  LQR  performance,  the  f^ure  of  the  NLQR  in 
this  situation  is  catastrophic,  leading  to  extremely  high  bandwidth  commands  and  raising 
the  likelihood  of  numerical  instabilities.  As  an  example.  Figure  5-11  contains  a  NLQR/EKF 
design  on  the  verge  of  instability.  As  such,  one  has  no  choice  but  to  increase  the  design 
envelope  to  encompass  all  expected  and  tolerable  velocities  reported  by  the  filter.  Toward 
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this  end,  the  limiting  of  ramp  velocity  during  maneuvers  becomes  not  just  a  good  rule  of 
thumb;  for  the  NLQR,  it  becomes  imperative. 

Because  the  NLQR  is  by  definition  accurate  over  a  larger  operating  regime  than  the 
LQR,  the  disadvantages  of  a  large  envelope  are  less  serious  for  the  NLQR.  The  primary 
difficulty  with  the  LQR  is  that  performance  must  always  be  traded  with  limit-cycling.  In 
Figure  5-15,  for  example,  the  loopshaping  response  limit  cycle  corresponds  to  fifteen-meter 
peak  to  peak  oscillations  for  the  full-scale  system.  The  designs  with  nonlinear  filtering 
and  the  NLQR  are  seen  to  have  less  limit  cycling,  especially  with  the  EKF.  In  fact,  direct 
comparisons  of  Figure  5-15  with  Figure  5-9,  and  of  Figure  5-16  with  Figure  5-10  show  that 
the  nonlinear  designs  with  an  envelope  of  0.4  meters  per  second,  are  superior  overall  to  the 
linear  designs  with  an  envelope  of  0.1  meters  per  second. 

Aside  from  the  limit  cycling  present  in  large-envelope  linear  designs,  the  linear  ap¬ 
proaches  prove  entirely  inappropriate  when  the  ramp  speed  is  increased.  For  instance. 
Figure  5-17  gives  only  the  responses  of  the  LQR/EKF,  NLQR/CGEKF,  and  NLQR/EKF 
controllers  because  the  other  designs  led  to  commands  that  were  beyond  the  capabilities  of 
the  apparatus.  Several  times,  the  cart  hit  the  side  of  the  tank,  and  in  all  cases,  cart  speeds 
of  greater  than  one  meter  per  second  were  demanded  by  the  controller.  Had  the  apparatus 
been  able  to  support  these  motions,  the  response  may  well  have  been  found  to  be  stable. 
Nonetheless,  the  ramp  speed  is  reasonable  at  0.4  meters  per  second,  so  the  linear  designs  are 
too  aggressive  at  best.  In  contrast,  the  motions  in  Figure  5-17  indicate  that  the  controllers 
are  able  to  execute  effective  overshoots,  reminiscent  of  our  results  in  Chapter  3. 

To  summarize,  the  use  of  purely  linear  designs  may  be  unsatisfactory  because  eliminat¬ 
ing  limit-cycling  sacrifices  performance,  and  improving  transient  performance  causes  limit 
cycling.  The  various  nonlinear  design  techniques  are  attractive  in  that  a  higher  design 
envelope  can  be  used  to  ensure  good  transient  performance,  but  with  reduced  oscillations 
because  no  linearizations  are  made.  Among  these  techniques,  the  LQR/EKF  appears  to 
be  the  best,  as  the  design  for  a  0.4  meters  per  second  envelope  worked  well  for  all  the 
runs.  The  NLQR/EKF  designs  were  equally  effective,  but  apparently  offer  no  significant 
improvements  over  the  LQR/EKF,  based  on  these  data.  Computationally,  the  NLQR  re¬ 
quires  the  construction  of  the  auxiliary  state  vector  at  each  time  step,  which  is  relatively 
insignificant  compared  to  propagation  of  the  EKF. 
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Figure  5-18:  Side  view  of  typical  step  responses. 

5.4.1  Heave  Response 

The  heave  response  of  the  plant  during  these  runs  is  not  large,  as  is  evident  in  Figure  5- 
18.  The  data  here  are  from  the  r\ms  of  Figure  5-14.  Vigorous  motions  at  the  top  of  the 
cable  agadn  lead  to  strong  vertical  excursions  at  the  towfrsh.  In  this  particular  case,  the 
linear  feedback  leads  to  vertical  excursions  that  are  three  times  as  large  as  for  the  nonlinear 
controllers.  As  stated  earlier,  these  excursions  should  not  be  taken  as  a  direct  scaling  of  the 
vertical  response  of  a  full-scale  system;  the  plot  is  intended  only  as  an  general  indication  of 
the  system’s  vertical  response. 
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5.4.2  Effects  of  Ambient  Currents  on  Oscillatory  Behavior 

The  developments  of  Chapter  2  indicate  that  the  lateral  frequency  response  of  a  cable  in 
water  is  strongly  attenuated  in  the  presence  of  ambient  currents.  Furthermore,  we  expect 
closed-loop  limit  cycling  to  be  reduced  in  currents,  based  on  the  (less  conclusive)  perturba¬ 
tion  results  of  that  chapter.  Here,  experimental  verification  of  this  hypothesis  is  considered, 
by  executing  feedback  dming  a  very  slow  setpoint  trajectory.  The  design  of  Figure  5-13  was 
incorporated  into  two-meter  moves  with  a  ramp  speed  of  0.05  meters  per  second,  and  the 
data  are  presented  in  Figure  5-19,  for  the  loopshaping  and  NLQR/EKF  designs.  In  these 
runs,  the  feedback  control  was  enabled  for  thirty  seconds  before  the  setpoint  move,  in  order 
to  allow  any  limit  cycles  to  develop  fully. 

Clearly,  the  ramp  speed  is  slow  enough  that  the  cable  in  both  cases  reaches  a  steady 
towing  configuration.  For  the  loopshaping  design,  the  nominal  limit  cycle  of  0.25  meters 
(peak-to-peak)  is  entirely  eliminated  during  the  transit,  but  then  reemerges  when  the  arti¬ 
ficial  current  is  “turned  off.”  The  NLQR/EKF  response  shows  no  limit  cycling  tendencies 
during  any  portion  of  the  nm,  as  expected. 

An  important  note  is  that  this  ramp  speed  is  very  small  on  an  ocean  scale,  since  the 
full-scale  speed  of  eight  centimeters  per  second  is  common  in  many  large  gyres  and  circula¬ 
tion  systems.  As  a  result,  many  full-scale  deployments  automatically  have  this  stabilizing 
element,  which  qualifies  the  linear  designs  to  some  extent. 

5.4.3  Incorporation  of  Actuator  Dynamics 

In  a  very  deep  deployment,  the  actuation  bandwidth  is  usually  very  much  higher  than 
that  of  the  cable  in  the  water.  In  the  present  case,  stiction  in  the  cart  drive  mechanism 
led  to  an  approximate  piue  delay  of  540  milliseconds,  which  corresponds  to  about  eleven 
seconds  in  the  full-scale  system.  These  actuation  dynamics  are  not  troublesome,  but  it 
may  make  sense  to  account  for  them  in  situations  where  the  delay  is  longer,  or  the  cable 
response  is  faster.  Toward  this  end,  one  could  append  the  plant  model  with  one  or  two 
states  to  model  the  actuator,  although  the  addition  of  these  states  might  make  the  NLQR 
computationally  difficult.  For  example,  the  fifth-order  nonlinear  control  design  of  Chapter  4 
requires  the  storage  of  a  15,625-element  matrix,  while  the  seventh-order  design  involves  a 
117,659-element  matrix. 
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Figure  5-19:  Experimental  verification  of  limit  cycle  attenuation  in  ambient  current. 
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Figure  5-20:  Generalized  Smith  feedback  block  diagram. 

Another  approach  that  could  be  used  is  ddayed  Smith  feedback  around  the  compensator 
[93].  The  Smith  predictor  has  been  used  extensively  in  linear  systems,  and  may  be  derived 
for  the  nonlinear  case  by  using  the  Loop  Transformation  Theorem  [33].  A  schematic  view 
of  generic  Smith  feedback  is  given  in  Figure  5-20.  The  basic  idea  is  first  to  construct 
the  compensator  K  to  stabilize  the  nominal  plant  Gnt  and  then  feed  back  around  the 
compensator  the  difference  between  the  outputs  of  the  nominal  plant  and  the  actual  plant 
model  Ga.  In  the  case  where  Ga  is  simply  Gn  with  a  pure  delay,  we  feed  back  a  delayed 
and  nondelayed  version  of  the  plant  output.  As  such,  from  a  computational  point  of  view, 
the  predictor  may  be  simply  added  onto  an  existing  compensator,  but  does  require  its  own 
state  estimation.  The  estimates  inherent  in  the  CGEKF  and  EKF  cannot  be  used  because 
they  are  continually  being  corrected. 

The  present  experiments  incorporate  the  delayed  Smith  scheme,  largely  because  the 
motor  servo  looks  more  like  a  pure  delay  (stiction  driven  by  the  integrator  in  the  motor 
.ervo)  than  a  first-  or  second-order  lag.  In  addition,  the  acoustic  fixes  arrive  one  time  step 
late,  so  setting  the  delay  at  820  milliseconds  accounts  for  both  delays.  An  illustration  of 
the  effects  of  delayed  Smith  feedback  is  given  in  Figure  5-21,  for  the  controllers  of  Figure  5- 
13.  The  transient  responses  appear  to  be  improved  only  slightly  for  these  runs,  while  the 
limit  cycle  is  practically  unchanged.  The  actuation  with  Smith  feedback  is  somewhat  less 
vigorous,  however,  and  a  larger  delay  would  presumably  enhance  the  stabilizing  properties 
of  this  scheme. 

The  remainder  of  tests  in  the  thesis  are  made  using  delayed  Smith  feedback,  including 
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Figure  5-21:  Effects  of  delayed  Smith  feedback  (0.82  seconds)  around  the  compensator. 

those  of  Figures  5*9  to  5-17  above.  An  important  note  is  that  for  linear  controllers,  the 
Smith  feedback  uses  a  linear  model.  In  addition,  when  errors  are  purposely  inserted  in  the 
nominal  plant  model,  these  errors  also  are  incorporated  in  the  Smith  feedback. 

5.4.4  Full  Smith  Feedback 

Another  situation  in  which  the  concept  of  Smith  feedback  is  useful  is  in  the  NLQR  designs 
with  a  large  envelope.  As  seen  in  Section  5.3.1,  the  cubic  approximation  of  quadratic  drag 
is  linear  near  the  origin;  on  the  tmit  line,  the  slope  at  the  origin  is  approximately  0.3. 
As  a  result,  when  the  envelope  is  large,  the  NLQR  may  induce  limit  cycling  for  the  same 
reasons  that  the  linear  designs  do.  Smith’s  scheme  can  be  used  effectively  in  this  case  also, 
by  letting  Ga  be  a  quadratic  model  with  delay,  and  letting  Gn  be  a  model  with  the  cubic 
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Figiire  5-22:  Effects  of  full  Smith  feedback  on  limit-cycling  behavior  of  the  NLQR/CGEKF. 

approximation  and  no  delay. 

This  technique  has  very  little  theoretical  justification,  since  the  model  with  quadratic 
drag  is,  of  course,  in  error.  Nonetheless,  it  is  presumably  closer  than  the  cubic  model. 
An  heuristic  attempt  was  made  to  reduce  the  oscillations  of  the  NLQR/CGEKF  design 
in  Figure  5-15,  which  has  an  envelope  of  0.4  meters  per  second,  and  oscillations  of  thirty 
centimeters  peak-to-peak  (ten  meters  in  the  full-scale  system).  Results  without  and  with 
the  full  Smith  feedback  are  given  in  Figure  5-22.  While  holding  station,  the  limit  cycle 
has  been  reduced  to  eighteen  centimeters  (six  meters  in  the  full-scale  system)  peak  to-peak. 
Furthermore,  we  see  that  the  oscillations  during  the  transit  are  eliminated  under  the  full 
Smith  feedback.  Finally,  an  artifact  of  this  strategy  has  the  towhsh  leading  the  setpoint, 
contrary  to  all  of  the  other  runs  made  in  this  chapter. 
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5.4.5  Disturbance  Rej^^on 

The  rejection  of  slow  disturbances  by  the  closed-loop  system  is  important  because  of  currents 
that  may  act  along  the  cable.  Simulating  this  situation  in  the  test  tank  proved  to  be 
problematic,  since  the  small  diameter  of  the  tank  made  it  impossible  to  reach  steady-state 
at  any  realistic  velocity.  A  simple  alternative  scheme  was  devised,  involving  a  nonlinear 
spring  attached  to  the  cable  midpoint:  see  Figure  5-23.  The  spring  consisted  of  a  weight 
and  a  string,  attached  to  the  tank  and  to  the  cable,  so  that  moving  away  from  the  walls 
caused  a  large  offset  in  the  cable  position.  Figure  5-24  shows  that  the  controllers  handle 
the  disturbance  as  expected.  The  design  was  the  same  as  in  Figure  5-13. 

Other  disturbances  which  are  likely  to  be  encormtered  in  ocean  operations  are  those  of 
wave  forces  on  the  surface  vessel,  and  limit  cycling  of  the  vessel  due  to  its  own  closed-loop 
dynamics.  Clearly,  the  wave-induced  motions  are  much  faster  than  the  limit  cycle,  and  they 
are  both  much  faster  than  the  cable.  In  order  to  assess  the  possible  impacts  on  stability 
and  performance  of  such  motions,  a  sinusoidal  cart  motion  was  imposed  upon  the  step 
responses  of  Figure  5-13.  This  motion  corresponds  to  a  ten-meter  peak-to-peak  surge  at 
.15  radians  per  second;  physical  constraints  prevented  the  scale  implementation  of  larger  or 
faster  motions.  Two  results  are  given  in  Figure  5-25;  clearly  there  is  no  loss  of  performance  or 
stability  associated  with  these  high-frequency  motions.  In  fact,  for  the  loopshaping  design, 
the  high-frequency  motions  eliminated  the  limit  cycling  almost  entirely,  and  they  reduced 
the  oscillations  in  the  NLQR/EKF  application.  This  trend  is  indicative  of  increased  system 
damping  due  to  high  velocities  in  the  upper  portion  of  the  cable.  For  both  controllers,  the 
motion  which  reached  the  bottom  of  the  cable  corresponds  to  approximately  1.5  meters 
peak-to-peak  in  the  full-scale  system.  Transmitted  motions  are  attenuated  further  when 
the  amplitudes  are  smaller,  or  when  the  frequency  is  higher. 

5.4.6  Higher-Order  Controller  Designs 

The  four-node  model  of  Figure  5-2  is  more  accurate  than  the  two-node  model,  for  the 
scaled  drill  string  system.  More  generally,  one  would  expect  that  higher-order  models  have 
even  better  performance  for  very  long  cables,  in  which  higher  modes  may  be  excited  and 
sustained.  Thus,  it  makes  sense  to  compare  the  performance  of  controllers  based  on  the 
four-node  model  with  those  of  the  two-node  model. 
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Figure  5-23:  Physical  application  of  an  unmodeled  spring  force  to  the  cable 
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Figure  5-24:  Illustration  of  slow  disturbance  rejection. 
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Figure  5-25:  Closed-loop  responses  in  the  presence  of  fast  cart  surge  motions. 
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The  design  programs  are  written  to  be  transparent  to  model  order,  so  the  only  difficulty 
is  the  construction  of  the  NLQR,  which  entails  a  model  order  of  length  nine  (including 
one  integrator),  a  lexicographic  state  vector  of  length  165,  and  a  triple  Kronecker  product 
with  531,441  elements.  With  double  precision  arithmetic,  this  large  matrix  temporarily 
occupies  approximately  4.25  megabytes  of  memory.  In  the  execution  of  the  program,  we 
simultaneously  propagated  a  ninth-order  EKF,  a  ninth-order  Smith  feedback  estimator,  and 
applied  a  9x165  feedback  matrix.  The  necessary  integrations  were  made  with  a  second- 
order  Adams-Bashforth  scheme  at  ten  cycles  per  acoustic  fix,  and  the  discrete-continuous 
forms  of  the  CGEKF  and  EKF  were  utilized  throughout.  In  the  experiment,  the  control 
loop  rate  was  kept  at  0.275  seconds.  These  runs  represent  the  computational  extremes  for 
this  chapter. 

The  runs  for  the  controllers  based  on  a  four-node  model  are  given  in  Figure  5-26;  the 
designs  were  based  on  the  same  envelope  (0.2  meters  per  second)  and  ramp  speed  (0.2  meters 
per  second)  as  in  Figure  5-13.  The  effects  of  model  order  are  clearer  in  Figures  5-27  and 
5-28,  in  which  the  tracking  error  of  the  towfish  is  plotted  for  the  the  two-node  and  four-node 
design  runs,  respectively.  The  integral  square  error  for  each  response  is  given  in  the  plots  as 
well.  We  see  that  the  four-node  loopshaping  result  is  very  poor,  while  the  LQR/KF  results 
are  virtually  identical  in  all  respects.  The  transient  responses  of  the  nonlinear  designs  in 
general  are  slightly  slower  for  the  four-node  designs,  but  have  fewer  residual  oscillations, 
and  greatly  reduced  limit  cycle  amplitudes.  The  reduction  of  oscillatory  behavior  in  the 
four-node  design  results  is  consistent  with  the  notion  that  higher-order  controller  designs 
can  account  for  more  structural  modes.  In  the  present  case,  our  familiar  limit  cycle  is 
eliminated  by  the  increased  resolution  of  the  four-node  controllers. 

5.4.7  Robustness 

In  this  section,  we  experimentally  investigate  the  robustness  properties  of  the  various  con¬ 
troller  design  techniques.  There  are  many  sources  of  model  uncertainty  in  a  towed  cable 
system  which  make  robustness  of  the  closed-loop  system  crucial.  To  begin,  the  effective  drag 
coefficient  along  the  cable  may  vary  a  great  deal,  due  to  complex  hydrodynamic  processes. 
For  example,  the  onset  of  vortex-induced  vibration  can  cause  the  cable  to  vibrate  in  the 
out-of-plane  direction  enough  to  double  the  effective  drag  coefficient.  These  vibrations  may 
be  steady  or  temporary,  and  may  be  localized  or  widespread,  or  even  travel  along  the  cable 
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Figure  5-26:  Step  responses  for  controllers  based  on  a  four-node  plant  model. 
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ire  5-27:  Tracking  errors  for  controllers  based  on  a  two-node  plant  model. 


Figure  5-28:  Tracking  errors  for  controllers  b£ised  on  a  four-node  plant  model. 
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[49].  At  the  other  extreme,  above  Re  ^  10^  one  may  encounter  supercritical  flow,  which 
leads  to  reduction  of  the  drag  coefficient  by  up  to  seventy  percent  [81j. 

Other  sources  of  modeling  error  could  include  possible  changes  in  the  unit  mass  of  a 
riser,  depending  on  what  materials  are  inside  (e.g.,  mud,  water,  slurry).  The  tensions  in 
such  a  riser  can  be  affected  if  bottom  contact  is  made,  or  if  the  payload  is  altered.  The 
tensions  are  also  changed  if  the  length  of  the  towed  system  is  modified;  ideally,  we  want  a 
controller  to  work  for  a  range  of  depths  near  the  design  depth. 

Linearization  and  the  Nature  of  Quadratic  Drag 

The  nonlinear  optimal  controller  of  Chapter  4  is  to  provide  a  control  law  based  on  nonlinear 
drag,  and  therefore  is  nominally  better  than  the  linear  control  law  for  a  range  of  the  state 
space.  A  velocity  envelope  is  selected  for  the  NLQR,  over  which  the  polynomial  approxi¬ 
mation  is  to  be  good,  and  the  linearization  for  linear  designs  is  selected  to  be  good  over  the 
same  domain.^  See,  for  example.  Figure  5-29. 

We  first  consider  the  case  in  which  the  quadratic  drag  coefficient  is  estimated  too  low; 
the  cubic  approximation  is  superior  to  the  linearization  near  the  origin  and  at  high  velocities 
both  inside  and  beyond  the  envelope.  See  Figure  5-29.  Figure  5-30  shows  a  case  where  the 
modeled  drag  coefficient  has  been  made  three  times  too  small.  The  design  parameters  for 
this  run  are  the  same  as  in  Figure  5-13;  the  bandwidth  is  0.20  radians  per  second,  with  an 
envelope  of  0.2  meters  per  second.  This  type  of  model  error  leads  to  very  stable  responses 
all  around,  as  the  low  product  of  linejirization  speed  and  modeled  drag  eliminates  limit 
cycling  in  the  linear  designs.  It  is  apparent  that  the  transient  performance  of  all  the  designs 
is  deteriorated  by  the  modeling  error,  although  less  for  the  nonlinear  designs  than  for  the 
linear  ones. 

In  the  event  that  the  drag  coefficient  is  overestimated,  the  properties  of  the  cubic  and 
linear  approximations  are  quite  different.  In  particular,  the  linearization  becomes  better 
than  the  cubic  approximation  for  velocities  in  the  higher  portion  of  the  domain,  and  for 
all  velocities  beyond  the  domain.  Thus,  in  this  case,  the  linear  controller  designs  may 
be  superior  because  the  nominal  linearization  is  closer  to  the  perturbed  quadratic  drag 
than  is  the  nominal  cubic  approximation.  This  fact  clearly  does  not  detract  from  the 

^Changing  the  linearization  speed  and  modifying  the  modeled  drag  coefficient  are  identical  from  the  point 
of  view  of  the  linear  designs,  but  quite  distinct  for  the  nonlinear  designs. 
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Figure  5*29:  Approximation  of  quadratic  drag  by  linearization  and  a  cubic  polynomial. 
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Figme  5*30:  Effects  of  underestimation  of  drag  coefficient. 
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Figure  5-31:  Effects  of  threefold  overestimation  of  drag  coefficient. 

fundamental  robustness  properties  of  the  nonlinear  optimal  regulator,  but  instead  enhances 
the  robustness  properties  of  the  linear  designs.  Safonov's  criteria  suggest  no  guarantees  for 
the  LQR  or  the  CGEKF  when  the  linearizations  are  too  high. 

In  the  experiments,  we  indeed  find  that  the  NLQR  caimot  provide  better  robustness  to 
this  type  of  drag  error  than  the  linear  designs;  see  Figme  5-31.  Here,  the  drag  coefficient 
was  overestimated  by  a  factor  of  three  in  the  design  plant  model,  the  envelope  was  0.2 
meters  per  second,  and  the  crossover  frequency  was  0.15  radians  per  second.  None  of  these 
designs  have  acceptable  performance,  as  the  responses  are  dominated  by  very  large  limit 
cycling.  The  corresponding  cart  motions  occupied  the  entire  tank,  and  the  cart  hit  the  wall 
fifty-five  seconds  into  the  LQR/KF  run. 

In  the  face  of  these  disappointing  responses,  there  are  several  redeeming  points  that  are 
noteworthy.  First,  the  controller  design  would  rarely  use  such  a  poor  model;  if  identification 
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runs  are  to  be  made,  these  would  ideally  incorporate  a  regime  of  supercritical  flow.  In 
addition,  the  dramatic  changes  in  drag  coefficient  occur  during  fully-developed  supercritical 
flow.  As  such,  these  changes  are  unlikely  to  he  seen  in  transient  motions  or  in  limit-cycling, 
although  an  exception  could  arise  during  operations  in  extremely  fast,  steady  currents 
(greater  than  one  meter  per  second). 

Another  point  is  that  the  onset  of  supercritical  flow  is  strongly  influenced  hy  the  rough¬ 
ness  of  the  cylinder.  For  example,  when  the  characteristic  roughness  ratio  (grain  size  divided 
hy  diameter)  is  0.005,  the  reduction  in  drag  coefficient  is  only  about  forty  percent  [81].  Sim¬ 
ilarly,  the  presence  of  kill-  and  choke-lines  on  a  marine  riser  could  be  advantageous.  Finally, 
we  reiterate  that  as  the  natural  plant  is  very  stable,  one  can  always  find  a  stabilizing  control 
law  by  reducing  the  bandwidth  or  speed  envelope  in  the  design. 

Robustness  to  Tension  and  Mass  Variations 

We  now  investigate  the  robustness  of  the  various  designs  with  specific  regard  to  tension  and 
mass  variations.  For  the  first  example,  the  weight  in  water  of  the  towfish  was  doubled,  as  was 
the  weight  of  the  cable,  in  the  design  plant  model.  As  before,  designs  with  a  speed  envelope 
of  0.2  meters  per  second  were  considered,  having  a  bandwidth  of  0.20  radians  per  second. 
The  results  are  shown  in  Figure  5-32.  The  fact  that  all  of  the  responses  seem  quite  stable 
is  not  surprising,  as  the  actual  plant  has  a  significantly  higher  natural  frequency  than  the 
control  law  expects.  In  comparison  with  Figure  5-14,  we  see  that  the  limit  cycling  is  almost 
completely  eliminated,  and  that  all  the  designs  suffer  some  loss  of  transient  performance. 
The  tracking  errors  are  plotted  for  the  two  cases  in  Figure  5-33  and  5-34,  for  the  the  correct 
and  incorrect  designs  respectively. 

Modeling  the  plant  as  too  light  was  expected  to  have  a  detrimental  effect  on  the  closed- 
loop  responses,  as  in  this  case  the  controller  is  paired  with  a  plant  whose  natural  frequencies 
are  reduced.  To  investigate  this  hypothesis,  the  design  of  Figure  5-13  was  incorporated,  but 
the  foam  flotation  on  the  cable  was  removed,  and  a  2.0-Newton  (in  water)  weight  was 
affixed  approximately  35  centimeters  up  from  the  tremsceiver.  Thus,  the  static  tensions 
were  increased  by  2.52  Newtons  uniformly.  The  results  are  shown  in  Figxire  5-35,  which 
may  be  compared  to  Figure  5-13,  generated  with  the  correct  model.  Overall,  this  modeling 
error  induces  a  limit  cycle  in  the  responses  of  the  nonlinear  designs,  but  does  not  affect  the 
11062^  responses  nearly  as  much.  The  amplitude  of  the  limit  cycle  in  the  NLQR/EKF  nm, 
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Figure  5-32;  Effects  of  designing  with  a  heavy  design  plant  model. 
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Figiire  5'33:  Tracking  errors  with  correct  design  plant  model. 


Figure  5-34:  Tracking  errors  with  heavy  design  plant  model. 


Figvire  5-35:  Effects  of  adding  weight  to  the  bottom  of  the  cable. 

for  example,  is  thirty-six  centimeters  peak-to-peak,  corresponding  to  approximately  twelve 
meters  for  the  full-scale  system.  There  is  no  associated  limit  cycle  in  5-13. 

As  all  the  designs  appear  susceptible  to  this  behavior,  our  only  consolation  is  that  these 
data  were  obtained  for  a  case  in  which  the  unmodeled  weight  was  equal  to  the  weight  of 
the  entire  modeled  system.  Obviously,  this  is  a  very  severe  error. 


5.5  Setpoint-Following 

The  good  results  of  the  previous  section,  and  those  of  Chapter  3,  suggest  the  design  of 
a  controller  with  the  necessary  capability  of  regulation  and  the  benefits  of  preshaping. 
More  specifically,  we  assume  that  a  vessel  trajectory  has  been  found  through  our  inversion 
technique,  and  now  we  would  like  to  apply  feedback  to  correct  for  modeling  errors  and 
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disturbances.  It  is  apparent  that  a  regulator  can  be  applied  to  the  plant  in  a  reference 
system  which  moves  along  the  preshaped  path.  In  other  words,  the  position  error  used  by 
the  regulator  is  the  distance  between  the  measured  towfish  and  the  moving  setpoint.  At  the 
same  time,  the  control  action  coming  out  of  the  regulator  is  added  to  the  moving  reference 
trajectory,  which  is  the  result  of  the  inversion. 

Our  first  example  is  given  for  the  scale  ARGO/JASON  system,  with  a  linearization 
speed  of  0.1  meters  per  second,  a  crossover  frequency  of  0.3  radians  per  second.  The  best 
two-node  model  available  was  used,  and  the  Smith  feedback  scheme  above  was  fully  utilized. 
For  the  preshaping,  a  towfish  trajectory  was  defined  which  involved  very  large  transients, 
as  well  as  periods  of  zero  and  low  velocity.  Figure  5-36  gives  the  results  for  five  different 
approaches  to  the  problem.  The  top  plot  shows  the  actual  motions  of  the  system,  while  the 
bottom  plot  indicates  the  tracking  errors.  The  top  curves  in  the  upper  plot  show  the  pure 
feedforward  solution,  with  no  feedback.  In  the  two  middle  sets  of  curves,  we  have  attempted 
to  follow  the  setpoint  using  the  loopshaping  and  NLQR/EKF  regulators.  Finally,  in  the 
bottom  pair  of  curves,  the  linear  and  nonlinear  regulators  have  been  layered  on  top  of  the 
preshaped  trajectories. 

In  general,  we  see  that  the  dynamic  inversion  for  this  run  is  very  successful.  In  fact,  it 
is  so  good  that  the  application  of  feedback  makes  the  response  slightly  worse.  As  expected, 
the  pure  regulation  results  are  qiiite  poor  due  to  the  conservativeness  in  the  designs.  The 
basic  conclusion  can  be  stated  as  follows:  when  the  inversion  is  known  to  he  good,  as  in 
the  cases  when  the  model  is  very  good  or  there  are  no  disturbances,  one  cannot  obtain 
better  performance  than  by  application  of  the  preshaped  input  alone.  This  conclusion  is 
not  surprising,  since  the  control  of  completely  known  plants  obviously  does  not  require 
feedback. 

We  now  consider  the  same  trajectory,  and  the  same  regulators  in  runs  where  the  non¬ 
linear  spring  of  Section  5.4.5  is  employed.  The  data  are  shown  in  Figure  5-37.  In  the 
application  of  the  preshaped  input  alone,  the  towfish  trajectory  first  shows  the  effects  of 
the  spring  (rather  suddenly)  at  approximately  eighteen  seconds  elapsed.  For  the  middle 
portion  of  the  run,  it  then  assumes  a  somewhat  constant  offset  from  the  setpoint,  crossing 
over  it  to  take  a  bias  in  the  other  direction  when  the  setpoint  goes  negative.  The  pure  reg¬ 
ulation  results  are  as  slow  as  in  the  imperturbed  case,  but  the  advantage  of  the  integrator 
in  the  control  is  evident  as  the  tracking  error  is  small  in  the  flat  middle  portion  of  the  run. 
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Figure  5-36:  Tracking  performance  in  the  absence  of  disturbances;  trajectories  (top)  and 
tracking  errors  (bottom). 


The  good  properties  of  both  the  preshaping  and  regulation  techniques  are  culminated  in  the 
final  two  curves,  where  the  control  is  layered  as  before.  The  feedback  becomes  active  only 
after  a  tracking  error  is  detected,  allowing  the  high  frequency  move  with  overshoot  to  be 
completed  successfully.  Then,  the  offset  due  to  the  disturbance  is  eliminated  by  feedback  on 
a  slower  time  scale.  The  layering  of  preshaping  trajectories  and  feedback  laws  is  thus  seen  to 
be  an  effective  means  for  achieving  the  high-bandwidth  trajectories  allowed  by  preshaping, 
and  the  important  corrective  action  of  feedback. 


5.6  Summary 

The  arguments  for  nonlinear  control  in  this  application  revolve  around  the  choice  of  design 
envelope,  and  the  resulting  dynamic  range.  The  data  show  that  the  linear  designs  can 
be  quite  stable  with  a  low  design  speed,  and  can  have  good  performance  with  a  high 
design  speed.  However,  it  seems  difficult  to  attain  both  properties  with  a  single  design,  for 
any  reasonable  range  of  operating  conditions.  The  nonlinear  designs,  in  contrast,  have  an 
improved  dynamic  range  which  makes  large-envelope  designs  successful  for  a  wide  operating 
regime.  This  improved  dynamic  range  is  especially  crucial  to  the  NLQR,  which  must  have 
a  large  envelope  anyway  to  enstire  that  the  polynomial  approximation  is  valid.  Overall, 
the  LQR/EKF  emerges  as  the  most  useful  among  the  four  nonlinear  approaches  that  were 
considered,  despite  the  predicted  failures  of  the  LQR  at  low  operating  speeds. 

From  a  design  point  of  view,  the  parameters  for  the  nonlinear  approaches  are  based 
wholly  on  those  of  linear  designs,  so  one  can  effectively  shape  the  nonlinear  controller  by 
manipulating  quantities  in  a  linear  framework.  Also,  the  auxiliary  scheme  of  delayed  Smith 
feedback  provides  a  usef\il  means  for  accommodating  pure  delays  in  the  control  loop,  while 
keeping  the  design  problem  tractable. 
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Figure  5-37:  Traddng  performance  in  the  presence  of  disturbances;  trajectories  (top)  and 
tracking  errors  (bottom). 


Chapter  6 


Conclusions  and 
Recommendations 

6.1  Summary  of  Chapter  2 

The  purpose  of  Chapter  2  was  to  develop  a  model  for  towed  cable  systems,  suitable  for 
simulation  and  controller  design,  and  to  examine  the  lateral  frequency  response  of  such  a 
cable.  Toward  this  end,  the  relevant  equations  describing  the  motions  of  underwater  cables 
were  described,  and  two  methods  of  approximate  analysis  were  presented.  The  harmonic 
balance  approach  gave  estimates  of  the  lateral  frequency  response  which  compared  well 
with  both  simulations  and  later  experimental  data.  This  response  is  typified  by  a  strong 
amplitude  dependence,  and  the  potential  loss  of  phase  at  arbitrarily  low  frequencies.  A 
multiple-scale  perturbation  analysis  was  also  conducted,  verifying  the  amplitude  dependence 
of  the  the  frequency  response,  and  providing  an  algebraic  relationship  between  the  excitation 
amplitude  and  frequency,  and  the  amplitude  of  the  towfish  motion.  This  compact  formula 
is  only  marginally  useful  away  from  the  natural  modes.  The  extension  of  this  perturbation 
approach  to  include  linear  feedback  was  successful  but  limited,  and  both  the  harmonic 
balance  and  perturbation  methods  correctly  predicted  the  attenuation  of  frequency  and 
closed-loop  responses  in  the  presence  of  currents  in  the  water. 
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6.2  Summary  of  Chapter  3 

The  ad  hoc  method  of  overshooting  a  desired  towiish  location  with  the  surface  vessel  was 
successfully  demonstrated  in  full-scale  tests,  with  a  potential  reduction  in  towiish  travel 
time  of  around  twenty  percent,  over  the  pure  step  responses.  A  cascade  of  out-of-plane 
overshoots  was  executed,  and  the  results  suggest  that  such  a  strategy  would  be  useful  for 
precision  grid-following  operations. 

The  idea  of  preshaping  vessel  trajectories  for  arbitrarily  complex  motions  of  the  towfish 
was  developed,  and  demonstrated  in  simulations  for  drill  strings  with  bending  stifbess, 
out-of-plane  towing  maneuvers,  and  the  heating  of  a  long  metal  ba^'  incorporating  nonlinear 
dynamics.  These  examples  indicate  the  generality  of  the  method,  and  illustrate  the  stability 
guidelines  for  the  technique.  The  method  was  successfully  demonstrated  in  full-scale  tests 
as  well,  which  punctuated  the  need  for  good  dynamic  positioning  of  the  vessel.  These 
tests  also  showed  also  that  since  there  is  no  feedback  involved  in  the  procedure,  navigation 
failures  cause  only  deterioration  of  the  towfish  trajectory.  This  makes  it  an  inherently  safe 
procedure,  which  is  very  attractive  from  an  operations  point  of  view. 


6.3  Summary  of  Chapter  4 

The  loop  operator  recovery  procedure  was  pursued  as  a  means  for  carrying  out  robust 
nonlinear  model-based  control.  As  this  procedure  benefits  from  nonlinear  observation  and 
nonlinear  feedback  laws,  a  number  of  estimators  were  outlined,  and  then  the  primary  con¬ 
tribution  of  the  chapter  was  presented:  a  practical  approximate  solution  to  the  nonlinear 
optimal  control  problem.  This  solution  is  valid  for  plant  models  having  analytic  nonlineeiri- 
ties  that  can  be  approximated  by  cubic  polynomials,  and  for  plants  having  arbitrary  order. 
Furthermore,  the  plant  may  be  multi-input /multi-output,  a  fact  which  should  make  the 
technique  useful  for  a  great  many  nonlinear  control  problems  beyond  this  thesis.  A  design 
example  for  a  sixth-order,  unstable  nonlinear  plant  illustrated  the  superior  performance  and 
robustness  properties  of  the  nonlinear  regulator,  compared  to  the  linear  quadratic  regulator. 


184 


6.4  Summary  of  Chapter  5 


Chapter  5  began  with  the  specification  of  two  full-scale  and  two  scale-model  systems,  and 
detailed  the  experimental  apparatus  at  the  test  tanlc  at  the  Woods  Hole  Oceanographic 
Institution.  One  setup  resembled  a  deep-water  ROV  system,  while  the  other  was  a  scaled 
drill  string  system.  Some  of  the  practical  issues  addressed  include  polynomial  approximation 
of  the  drag  nonlinearity,  quantification  of  modeling  errors  at  the  plant  input,  and  calculation 
of  the  guaranteed  properties  of  the  Constant  Extended  Kalman  Filter  and  the  linear 
quadratic  regulator. 

The  closed-loop  response  to  a  ramped  step  change  in  the  setpoint  was  adopted  as  the 
standard  test,  and  numerous  experiments  with  the  scaled  drill  string  system  indicated  that 
the  nonlinear  approaches  have  significant  advantages  over  linear  controllers.  In  particular, 
the  linear  approaches  appear  to  trade  limit-cycling  for  poor  performance,  depending  on  the 
linearization  chosen.  In  contrast,  the  nonlinear  approaches  have  a  much  larger  dynamic 
range,  so  that  one  can  achieve  good  performance  and  long-term  stability  simultaneously. 
This  d3mamic  range  is  fortunate,  as  it  was  illustrated  that  the  envelope  for  the  nonlinear 
optimal  control  law  must  encompass  the  allowed  operating  conditions,  or  else  catastrophic 
instability  can  occur. 

A  niimber  of  variations  were  investigated,  including  the  application  of  delayed  Smith 
feedback,  the  use  of  high-order  models  in  controller  design,  and  the  rejection  of  fast  and  slow 
disturbances  acting  on  the  system.  The  robustness  properties  of  the  designs  were  considered 
as  well,  and  the  linear  designs  were  found  to  be  remarkably  capable  in  this  regard.  More 
specifically,  when  stabilizing  errors  were  made  in  the  tension  and  drag  terms  of  the  design 
plant  model,  there  were  found  only  moderate  improvements  in  transient  performance  by 
using  the  nonlinear  approaches.  However,  when  destabilizing  errors  were  made,  the  failures 
of  all  the  designs  looked  very  similar.  This  trend  is  in  contradistinction  to  that  of  the 
design  example  in  Chapter  4,  where  the  linear  quadratic  reg^ator  performed  very  poorly. 
The  main  difference  between  that  example  and  the  later  cable  tests  is  that  in  the  latter 
case,  drag  is  naturally  stabilizing,  which  always  benefits  the  linear  designs. 

Finally,  the  layering  of  regulation  designs  onto  a  preshaped  reference  frame,  for  the  scaled 
ROV  system,  was  verified.  This  layered  control  approach  takes  advantage  of  the  preshaping 
during  high-bandwidth  motions,  and  of  the  feedback  during  slow  motions.  Therefore,  the 
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usual  rule  that  robustness  and  performance  cannot  be  achieved  simultaneously  appears 
to  be  validated.  Our  approach  provides  a  smooth  way  to  transition  between  regimes  of 
performance-oriented  operation,  and  robustness-oriented  operation. 


6.5  The  Big  Picture  and  Recommendations 

This  thesis  has  presented  techniques  that  are  useful  for  a  range  of  typical  deep-ocean  oper¬ 
ations  employing  long  vertical  cables.  Along  the  way,  more  general  procedures  for  dynamic 
inversion  and  for  nonlinear  optimal  regulator  design  were  developed.  The  inversion  scheme, 
as  it  stands,  has  direct  utility  for  many  towing  problems,  and  it  has  been  verified  repeatedly 
in  full-scale  and  model  tests,  and  in  simulations  for  various  other  nonlinear  systems.  The 
approach  shoiild  be  useful  for  a  range  of  tridiagonal  plants  with  distributed  states  and  local 
nonlinearities.  The  nonlinear  optimal  regvilator  solution  that  has  been  presented,  however, 
has  a  number  of  limitations.  To  begin,  as  was  stated  previously,  the  “hardening”  type  of 
nonlinearity  is  likely  to  be  most  appropriate  for  the  method,  as  the  core  of  the  procedure 
entails  making  polynomial  approximations  of  the  nonlinearities.  Furthermore,  it  was  found 
that  the  expansions  become  very  difficult  as  the  polynomial  order  exceeds  three.  Common 
saturation-  or  deadband-type  nonlineaiities  cannot  be  handled  easily  with  this  approach; 
perhaps  future  work  will  enable  the  use  of  different  basis  functions.  More  work  is  also 
needed  in  the  area  of  model-based  nonlinear  control,  as  the  LOR  technique  is  theoretically 
valid  only  when  the  controller  loop  is  recovered.  At  the  same  time,  the  experimental  results 
of  the  thesis  are  very  encouraging  for  nonlinear  control  of  this  type. 

Finally,  the  range  of  control  problems  associated  with  towed  imderwater  cables  is  much 
larger  than  that  covered  in  this  thesis.  Tests  need  to  be  run  in  the  presence  of  large  currents 
which  could  affect  the  closed-loop  properties  of  our  designs.  The  test  tank  was  inadequate 
for  these  purposes,  although  experiments  in  a  towing  tank  could  be  conducted.  In  addition, 
complications  can  arise  in  cables  with  large  layback,  due  to  geometric  coupling  between 
in-plane  and  out-of-plane,  as  well  as  vertical,  motions.  These  situations  will  require  models 
of  increased  complexity,  and  multi-input,  multi-output  control  schemes  of  commensurate 
ability. 
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Appendix  A 

Mathematical  Preliminaries 


In  nonlinear  control,  the  I/O  viewpoint  is  often  useful,  in  which  operators  map  signals  from 
an  input  space  to  an  output  space.  A  standard  reference  for  this  section  is  Desoer  and 
Vidyasagar  [33].  The  signals  are  entire  trajectories  and  not  just  points  in  R'*,  and  *heir  size 
is  usually  described  with  norms;  we  define  the  p-Dorm  of  a  signal  z  to  be 

11*11,=  [/“  1*1' <*1^  (A.1) 

In  the  case  that  p  =  oo,  we  have  ||  z  ||ao=  sup^  |  z(f)  |  .  Let  Lp  be  the  set  of  all  signals  for 
which  II  z  lip  is  finite,  i.e.. 


Lp={x\\\x  \\p<  oo} 


(A.2) 


The  truncation  operator  Pr  is  given  by 


(PrX)  =  { 


Z(t),  t  <  T 

0,  t  >  T 

and  the  extended  space  L‘  consists  of  all  signals  with  truncations  in  Lp,  i.e.. 


(A.3) 


X;  =  {z|P,z3Xp,Vr>0}.  (A.4) 

Thus,  X'  includes  signals  which  are  growing  without  bound,  as  well  as  those  which  are 
infinite  over  a  set  of  zero  measure.  Finally,  the  truncated  norm  of  a  signal  is  given  by 
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The  Lp-Dorm^  or  gain,  of  an  operator  g  is  defined  as 


(A.5) 


II  “  I|P.T 


(A.6) 


The  Lp  incremental  gain  of  g  is 


llflllp.<=  s"P 


II  ffUi  -  gM2  ||p.T 

ll«l-tt2llp,T  ■ 


(A.7) 


We  say  that  an  operator  is  Lp-stable  if  it  has  finite  gain;  this  property  is  often  referred  to  as 
finite-gain  stability.  Similarly,  if  inputs  of  bounded  norm  induce  outputs  of  bounded  norm, 
a  system  is  said  to  be  bounded. 

A  nonlinear  estimator  generates  an  estimate  x,  so  that  the  observation  error  e  :=  i  —  x 
can  be  viewed  as  the  output  of  an  operator  £!u((t  which  takes  as  its  inputs  the  process 
and  measurement  noises  and  depends  on  u.  The  estimator  is  said  to  be  nondivergent  [82] 
if  there  exists  a  contmuous  increasing  function  p:  R  such  that 


S«P  II  e  llp.T<  p{\\  9  ||p,r).  (A.8) 

^,-r 

More  informally,  a  nondivergent  estimator  takes  norm-bounded  disturbances  (  and  0  into 
norm-bounded  errors  in  the  estimate.  Not«  that  nondivergence  has  to  be  specified  with 
respect  to  a  given  norm.  The  estimate  is  nondivergent  with  finite  gain  if  the  norm-bound 
of  the  error  is  proportional  to  the  magnitude  of  the  disturbances. 
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